{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Back to the widget production problem\n",
    "Let's go back to the Opti 101/201 classic: Widget production and distribution. Use the code below as the base model for this problem."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%pip install gurobipy\n",
    "\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "import gurobipy as gp\n",
    "from gurobipy import GRB"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align*} \n",
    "{\\rm minimize} \\space &\\sum_{p,d}c_{p,d}x_{p,d}, \\forall p \\in P, d \\in D\\\\ \n",
    "\\sum_{p}&x_{p,d} \\ge n_d, \\forall d \\in D \\\\ \n",
    "\\sum_{d}&x_{p,d} \\le m_p, \\forall p \\in P \\\\ \n",
    "\\sum_{d}&x_{p,d} \\ge a*m_p, \\forall p \\in P\\\\ \n",
    "&x_{p,d} \\in \\{0\\} \\cup [30, u_p]\n",
    "\\end{align*}\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-warning\">\n",
    "    <strong>Note!</strong>\n",
    "    <p>Make sure to run the next two code cells to get a baseline for the original problem.</p>\n",
    "</div>\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# read in transportation cost data\n",
    "path = 'https://raw.githubusercontent.com/Gurobi/modeling-examples/master/optimization101/Modeling_Session_1/'\n",
    "transp_cost = pd.read_csv(path + 'cost.csv')\n",
    "\n",
    "# get production and distribution locations from data frame\n",
    "production = list(transp_cost['production'].unique())\n",
    "distribution = list(transp_cost['distribution'].unique())\n",
    "transp_cost = transp_cost.set_index(['production','distribution']).squeeze()\n",
    "\n",
    "max_prod = pd.Series([180,200,140,80,180], index = production, name = \"max_production\")\n",
    "n_demand = pd.Series([89,95,121,101,116,181], index = distribution, name = \"demand\")\n",
    "\n",
    "# min overall production requirement, and min shipment amount from p to d\n",
    "frac = 0.75\n",
    "C = 30"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Code for the 'gurobipy' model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "This model has a total cost of 1748.42\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th>cost</th>\n",
       "      <th>shipment</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>production</th>\n",
       "      <th>distribution</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th rowspan=\"2\" valign=\"top\">Baltimore</th>\n",
       "      <th>Nashville</th>\n",
       "      <td>5.96</td>\n",
       "      <td>30.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Richmond</th>\n",
       "      <td>1.96</td>\n",
       "      <td>116.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"2\" valign=\"top\">Cleveland</th>\n",
       "      <th>Columbia</th>\n",
       "      <td>2.43</td>\n",
       "      <td>89.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Indianapolis</th>\n",
       "      <td>2.37</td>\n",
       "      <td>95.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Little Rock</th>\n",
       "      <th>St. Louis</th>\n",
       "      <td>2.92</td>\n",
       "      <td>140.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Birmingham</th>\n",
       "      <th>Nashville</th>\n",
       "      <td>1.53</td>\n",
       "      <td>71.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"2\" valign=\"top\">Charleston</th>\n",
       "      <th>Lexington</th>\n",
       "      <td>1.61</td>\n",
       "      <td>121.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>St. Louis</th>\n",
       "      <td>4.60</td>\n",
       "      <td>41.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                          cost  shipment\n",
       "production  distribution                \n",
       "Baltimore   Nashville     5.96      30.0\n",
       "            Richmond      1.96     116.0\n",
       "Cleveland   Columbia      2.43      89.0\n",
       "            Indianapolis  2.37      95.0\n",
       "Little Rock St. Louis     2.92     140.0\n",
       "Birmingham  Nashville     1.53      71.0\n",
       "Charleston  Lexington     1.61     121.0\n",
       "            St. Louis     4.60      41.0"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m = gp.Model('widgets')\n",
    "m.setParam('OutputFlag',0)\n",
    "m.setParam('PoolSearchMode',2)   \n",
    "m.setParam('PoolSolutions',2)   \n",
    "\n",
    "# decision vars\n",
    "x = m.addVars(production, distribution, vtype=GRB.SEMICONT, lb = C, name = 'prod_ship')\n",
    "\n",
    "# constraints\n",
    "can_produce = m.addConstrs((gp.quicksum(x[p,d] for d in distribution) <= max_prod[p] for p in production), name = 'can_produce')\n",
    "must_produce = m.addConstrs((gp.quicksum(x[p,d] for d in distribution) >= frac*max_prod[p] for p in production), name = 'must_produce')\n",
    "meet_demand = m.addConstrs((x.sum('*', d) >= n_demand[d] for d in distribution), 'meet_demand')\n",
    "\n",
    "# objective\n",
    "total_cost = gp.quicksum(transp_cost[p,d]*x[p,d] for p in production for d in distribution)\n",
    "m.setObjective(total_cost, GRB.MINIMIZE)\n",
    "\n",
    "# optimize and extract and show solution\n",
    "m.optimize()\n",
    "\n",
    "x_values = pd.Series(m.getAttr('X', x), name = \"shipment\", index = transp_cost.index)\n",
    "sol = pd.concat([transp_cost, x_values], axis=1)\n",
    "\n",
    "m_ObjVal = m.ObjVal\n",
    "print(f\"This model has a total cost of {round(m_ObjVal,2)}\")\n",
    "sol[sol.shipment > 0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Piecewise Linear Modeling"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-success alert-dismissible\">\n",
    "  <a href=\"#\" class=\"close\" data-dismiss=\"alert\" aria-label=\"close\">&times;</a>\n",
    "    <strong>Here are a couple a good documentation links to help you use PWL in the code below. </strong>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Links to PWL in `gurobipy`:\n",
    "- [Model.addGenConstrPWL](https://docs.gurobi.com/projects/optimizer/en/current/reference/python/model.html#Model.addGenConstrPWL) command\n",
    "- [Visualization of PWL](https://docs.gurobi.com/projects/optimizer/en/current/concepts/modeling/objectives.html#piecewise-linear-objectives)\n",
    "\n",
    "### Scenario 1\n",
    "We have a new cost structure for transporting widgets we make. \n",
    "- There will be a flat 0.50 increase for all transportation costs. \n",
    "- If we decide to bulk ship at least half of a production facility's max to any one distribution location, every widget over that half capacity number costs 60% of the new transportation cost. \n",
    "\n",
    " So, if the original cost is 2, the new cost of shipping from $p$ to $d$ for the first half of capacity is 2.5 (0.50 flat increase), and anything produced after that ships at 1.5. Update the model to reflect this (0.60 * 2.5). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "This model has a total cost of 1925.29\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th>cost</th>\n",
       "      <th>base_costs</th>\n",
       "      <th>bulk_costs</th>\n",
       "      <th>shipment_pwl</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>production</th>\n",
       "      <th>distribution</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Baltimore</th>\n",
       "      <th>Richmond</th>\n",
       "      <td>1.96</td>\n",
       "      <td>2.46</td>\n",
       "      <td>1.48</td>\n",
       "      <td>135.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">Cleveland</th>\n",
       "      <th>Columbia</th>\n",
       "      <td>2.43</td>\n",
       "      <td>2.93</td>\n",
       "      <td>1.76</td>\n",
       "      <td>89.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Indianapolis</th>\n",
       "      <td>2.37</td>\n",
       "      <td>2.87</td>\n",
       "      <td>1.72</td>\n",
       "      <td>40.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Nashville</th>\n",
       "      <td>4.13</td>\n",
       "      <td>4.63</td>\n",
       "      <td>2.78</td>\n",
       "      <td>30.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>St. Louis</th>\n",
       "      <td>4.88</td>\n",
       "      <td>5.38</td>\n",
       "      <td>3.23</td>\n",
       "      <td>41.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Little Rock</th>\n",
       "      <th>St. Louis</th>\n",
       "      <td>2.92</td>\n",
       "      <td>3.42</td>\n",
       "      <td>2.05</td>\n",
       "      <td>140.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Birmingham</th>\n",
       "      <th>Nashville</th>\n",
       "      <td>1.53</td>\n",
       "      <td>2.03</td>\n",
       "      <td>1.22</td>\n",
       "      <td>71.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"2\" valign=\"top\">Charleston</th>\n",
       "      <th>Indianapolis</th>\n",
       "      <td>2.61</td>\n",
       "      <td>3.11</td>\n",
       "      <td>1.87</td>\n",
       "      <td>55.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Lexington</th>\n",
       "      <td>1.61</td>\n",
       "      <td>2.11</td>\n",
       "      <td>1.27</td>\n",
       "      <td>121.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                          cost  base_costs  bulk_costs  shipment_pwl\n",
       "production  distribution                                            \n",
       "Baltimore   Richmond      1.96        2.46        1.48         135.0\n",
       "Cleveland   Columbia      2.43        2.93        1.76          89.0\n",
       "            Indianapolis  2.37        2.87        1.72          40.0\n",
       "            Nashville     4.13        4.63        2.78          30.0\n",
       "            St. Louis     4.88        5.38        3.23          41.0\n",
       "Little Rock St. Louis     2.92        3.42        2.05         140.0\n",
       "Birmingham  Nashville     1.53        2.03        1.22          71.0\n",
       "Charleston  Indianapolis  2.61        3.11        1.87          55.0\n",
       "            Lexington     1.61        2.11        1.27         121.0"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m_pwl = gp.Model('widgets')\n",
    "m_pwl.setParam('OutputFlag', 0)\n",
    "\n",
    "x_pwl = m_pwl.addVars(production, distribution, vtype=GRB.SEMICONT, lb = 30, name = 'prod_ship')\n",
    "pwl_cost = m_pwl.addVars(production, distribution, name = 'pwl_cost')\n",
    "\n",
    "# constraints\n",
    "can_produce = m_pwl.addConstrs((gp.quicksum(x_pwl[p,d] for d in distribution) <= max_prod[p] for p in production), name = 'can_produce')\n",
    "must_produce = m_pwl.addConstrs((gp.quicksum(x_pwl[p,d] for d in distribution) >= frac*max_prod[p] for p in production), name = 'must_produce')\n",
    "meet_demand = m_pwl.addConstrs((x_pwl.sum('*', d) >= n_demand[d] for d in distribution), name = 'demand')\n",
    "\n",
    "# create empty series for analysis \n",
    "bulk_costs = pd.Series(index=transp_cost.keys(), name='bulk_costs')\n",
    "base_costs = pd.Series(index=transp_cost.keys(), name='base_costs')\n",
    "\n",
    "# Create piecewise cost functions\n",
    "for p in production:\n",
    "    for d in distribution:\n",
    "        max_produce = max_prod[p]\n",
    "        frac_produce = 0.5*max_produce\n",
    "        \n",
    "        # Adjusted costs\n",
    "        base_cost = transp_cost[p,d] + 0.5\n",
    "        bulk_cost = round(0.6*base_cost, 2) # 25% increase for the first half\n",
    "\n",
    "        base_costs[p,d] = base_cost\n",
    "        bulk_costs[p,d] = bulk_cost\n",
    "\n",
    "        # Define the piecewise-linear function points\n",
    "        pwl_points_x = [\n",
    "            0, frac_produce, max_produce\n",
    "        ]\n",
    "\n",
    "        pwl_points_y = [\n",
    "            0, base_cost * frac_produce, base_cost * frac_produce + bulk_cost * (max_produce - frac_produce)\n",
    "        ]\n",
    "\n",
    "        m_pwl.addGenConstrPWL(x_pwl[p, d], pwl_cost[p, d], pwl_points_x, pwl_points_y, name=f\"pwl_{p}_{d}\")\n",
    "\n",
    "# Objective: minimize the total cost\n",
    "m_pwl.setObjective(gp.quicksum(pwl_cost[p, d] for p, d in transp_cost.keys()), GRB.MINIMIZE)\n",
    "\n",
    "# Optimize\n",
    "m_pwl.optimize()\n",
    "\n",
    "x_values_pwl = pd.Series(m_pwl.getAttr('X', x_pwl), name = \"shipment_pwl\", index = transp_cost.index)\n",
    "sol_pwl = pd.concat([transp_cost, base_costs, bulk_costs, x_values_pwl], axis=1)\n",
    "\n",
    "print(f\"This model has a total cost of {round(m_pwl.ObjVal,2)}\")\n",
    "sol_pwl[sol_pwl.shipment_pwl > 0]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th>cost</th>\n",
       "      <th>shipment</th>\n",
       "      <th>base_costs</th>\n",
       "      <th>bulk_costs</th>\n",
       "      <th>shipment_pwl</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>production</th>\n",
       "      <th>distribution</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th rowspan=\"2\" valign=\"top\">Baltimore</th>\n",
       "      <th>Nashville</th>\n",
       "      <td>5.96</td>\n",
       "      <td>30.0</td>\n",
       "      <td>6.46</td>\n",
       "      <td>3.88</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Richmond</th>\n",
       "      <td>1.96</td>\n",
       "      <td>116.0</td>\n",
       "      <td>2.46</td>\n",
       "      <td>1.48</td>\n",
       "      <td>135.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">Cleveland</th>\n",
       "      <th>Columbia</th>\n",
       "      <td>2.43</td>\n",
       "      <td>89.0</td>\n",
       "      <td>2.93</td>\n",
       "      <td>1.76</td>\n",
       "      <td>89.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Indianapolis</th>\n",
       "      <td>2.37</td>\n",
       "      <td>95.0</td>\n",
       "      <td>2.87</td>\n",
       "      <td>1.72</td>\n",
       "      <td>40.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Nashville</th>\n",
       "      <td>4.13</td>\n",
       "      <td>0.0</td>\n",
       "      <td>4.63</td>\n",
       "      <td>2.78</td>\n",
       "      <td>30.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>St. Louis</th>\n",
       "      <td>4.88</td>\n",
       "      <td>0.0</td>\n",
       "      <td>5.38</td>\n",
       "      <td>3.23</td>\n",
       "      <td>41.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Little Rock</th>\n",
       "      <th>St. Louis</th>\n",
       "      <td>2.92</td>\n",
       "      <td>140.0</td>\n",
       "      <td>3.42</td>\n",
       "      <td>2.05</td>\n",
       "      <td>140.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Birmingham</th>\n",
       "      <th>Nashville</th>\n",
       "      <td>1.53</td>\n",
       "      <td>71.0</td>\n",
       "      <td>2.03</td>\n",
       "      <td>1.22</td>\n",
       "      <td>71.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"3\" valign=\"top\">Charleston</th>\n",
       "      <th>Indianapolis</th>\n",
       "      <td>2.61</td>\n",
       "      <td>0.0</td>\n",
       "      <td>3.11</td>\n",
       "      <td>1.87</td>\n",
       "      <td>55.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Lexington</th>\n",
       "      <td>1.61</td>\n",
       "      <td>121.0</td>\n",
       "      <td>2.11</td>\n",
       "      <td>1.27</td>\n",
       "      <td>121.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>St. Louis</th>\n",
       "      <td>4.60</td>\n",
       "      <td>41.0</td>\n",
       "      <td>5.10</td>\n",
       "      <td>3.06</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                          cost  shipment  base_costs  bulk_costs  shipment_pwl\n",
       "production  distribution                                                      \n",
       "Baltimore   Nashville     5.96      30.0        6.46        3.88           0.0\n",
       "            Richmond      1.96     116.0        2.46        1.48         135.0\n",
       "Cleveland   Columbia      2.43      89.0        2.93        1.76          89.0\n",
       "            Indianapolis  2.37      95.0        2.87        1.72          40.0\n",
       "            Nashville     4.13       0.0        4.63        2.78          30.0\n",
       "            St. Louis     4.88       0.0        5.38        3.23          41.0\n",
       "Little Rock St. Louis     2.92     140.0        3.42        2.05         140.0\n",
       "Birmingham  Nashville     1.53      71.0        2.03        1.22          71.0\n",
       "Charleston  Indianapolis  2.61       0.0        3.11        1.87          55.0\n",
       "            Lexington     1.61     121.0        2.11        1.27         121.0\n",
       "            St. Louis     4.60      41.0        5.10        3.06           0.0"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sol_comp = pd.concat([sol, sol_pwl.drop('cost', axis=1)], axis=1)\n",
    "sol_comp[(sol_comp.shipment > 0) | (sol_comp.shipment_pwl > 0)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Scenario 2\n",
    "The original model does not reflect any production costs to make a widget in different locations. Suppose we have a given production cost per widget for making up to 25% of a production facility's capacity. After that level of production is reached, then the cost *increases* by 25%. For example, if a facility can produce 100 widgets and the initial cost is $1, then the first 25 widgets made will cost $1 each and up to the next 75 widgets qill cost $1.25 each to make. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Gurobi Optimizer version 12.0.0 build v12.0.0rc1 (mac64[rosetta2] - Darwin 24.1.0 24B83)\n",
      "\n",
      "CPU model: Apple M1\n",
      "Thread count: 8 physical cores, 8 logical processors, using up to 8 threads\n",
      "\n",
      "Optimize a model with 21 rows, 40 columns and 125 nonzeros\n",
      "Model fingerprint: 0x4405fcbf\n",
      "Model has 5 simple general constraints\n",
      "  5 PWL\n",
      "Variable types: 10 continuous, 0 integer (0 binary)\n",
      "Semi-Variable types: 30 continuous, 0 integer\n",
      "Coefficient statistics:\n",
      "  Matrix range     [1e+00, 1e+00]\n",
      "  Objective range  [1e+00, 7e+00]\n",
      "  Bounds range     [3e+01, 3e+01]\n",
      "  RHS range        [6e+01, 2e+02]\n",
      "  PWLCon x range   [2e+01, 5e+01]\n",
      "  PWLCon y range   [6e+01, 2e+02]\n",
      "Presolve removed 10 rows and 5 columns\n",
      "Presolve time: 0.00s\n",
      "Presolved: 71 rows, 65 columns, 185 nonzeros\n",
      "Variable types: 35 continuous, 30 integer (30 binary)\n",
      "\n",
      "Root relaxation: objective 3.073585e+03, 27 iterations, 0.00 seconds (0.00 work units)\n",
      "\n",
      "    Nodes    |    Current Node    |     Objective Bounds      |     Work\n",
      " Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time\n",
      "\n",
      "     0     0 3073.58500    0    3          - 3073.58500      -     -    0s\n",
      "H    0     0                    3121.2950000 3073.58500  1.53%     -    0s\n",
      "     0     0     cutoff    0      3121.29500 3121.29500  0.00%     -    0s\n",
      "\n",
      "Cutting planes:\n",
      "  Gomory: 2\n",
      "  Cover: 1\n",
      "  Implied bound: 2\n",
      "  MIR: 3\n",
      "  Flow cover: 6\n",
      "  Relax-and-lift: 2\n",
      "\n",
      "Explored 1 nodes (44 simplex iterations) in 0.03 seconds (0.00 work units)\n",
      "Thread count was 8 (of 8 available processors)\n",
      "\n",
      "Solution count 1: 3121.29 \n",
      "\n",
      "Optimal solution found (tolerance 1.00e-04)\n",
      "Best objective 3.121295000000e+03, best bound 3.121295000000e+03, gap 0.0000%\n",
      "This model has a total cost of 3121.29\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th>shipment</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>production</th>\n",
       "      <th>distribution</th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th rowspan=\"2\" valign=\"top\">Baltimore</th>\n",
       "      <th>Nashville</th>\n",
       "      <td>30.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Richmond</th>\n",
       "      <td>116.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"2\" valign=\"top\">Cleveland</th>\n",
       "      <th>Columbia</th>\n",
       "      <td>89.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Indianapolis</th>\n",
       "      <td>95.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Little Rock</th>\n",
       "      <th>St. Louis</th>\n",
       "      <td>140.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Birmingham</th>\n",
       "      <th>Nashville</th>\n",
       "      <td>71.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"2\" valign=\"top\">Charleston</th>\n",
       "      <th>Lexington</th>\n",
       "      <td>121.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>St. Louis</th>\n",
       "      <td>41.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                          shipment\n",
       "production  distribution          \n",
       "Baltimore   Nashville         30.0\n",
       "            Richmond         116.0\n",
       "Cleveland   Columbia          89.0\n",
       "            Indianapolis      95.0\n",
       "Little Rock St. Louis        140.0\n",
       "Birmingham  Nashville         71.0\n",
       "Charleston  Lexington        121.0\n",
       "            St. Louis         41.0"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m_pwl2 = gp.Model('widgets')\n",
    "#m_pwl2.setParam('OutputFlag',0)\n",
    "\n",
    "# decision vars\n",
    "x_pwl2 = m_pwl2.addVars(production, distribution, vtype=GRB.SEMICONT, lb = C, name = 'prod_ship')\n",
    "\n",
    "# constraints\n",
    "can_produce = m_pwl2.addConstrs((gp.quicksum(x_pwl2[p,d] for d in distribution) <= max_prod[p] for p in production), name = 'can_produce')\n",
    "must_produce = m_pwl2.addConstrs((gp.quicksum(x_pwl2[p,d] for d in distribution) >= frac*max_prod[p] for p in production), name = 'must_produce')\n",
    "meet_demand = m_pwl2.addConstrs(x_pwl2.sum('*', d) >= n_demand[d] for d in distribution)\n",
    "\n",
    "# Total production for each facility\n",
    "base_prod_cost = pd.Series([1,1.2,.8,.8,.9], index = production, name = \"max_production\")\n",
    "\n",
    "# Add new variables and/or constraints\n",
    "q = m_pwl2.addVars(production, name='q')\n",
    "pwl_prod = m_pwl2.addVars(production, name='pwl_prod')\n",
    "m_pwl2.addConstrs((q[p] == gp.quicksum(x_pwl2[p, d] for d in distribution) for p in production), name='total_production')\n",
    "\n",
    "for p in production:\n",
    "    m_pwl2.addGenConstrPWL(\n",
    "        q[p],\n",
    "        pwl_prod[p],\n",
    "        [0, 0.25 * max_prod[p], max_prod[p]],\n",
    "        [0, base_prod_cost[p]*max_prod[p], base_prod_cost[p]*max_prod[p] + 1.25*base_prod_cost[p]* (max_prod[p] - 0.25*max_prod[p])],\n",
    "        name=f'pwl_cost_{p}'\n",
    "    )\n",
    "\n",
    "# objective\n",
    "total_transp_cost = gp.quicksum(transp_cost[p,d]*x_pwl2[p,d] for p in production for d in distribution)\n",
    "m_pwl2.setObjective(total_transp_cost + pwl_prod.sum(), GRB.MINIMIZE)\n",
    "\n",
    "# optimize and extract and show solution\n",
    "m_pwl2.optimize()\n",
    "\n",
    "x_values_prod = pd.Series(m_pwl2.getAttr('X', x_pwl2), name = \"shipment\", index = transp_cost.index)\n",
    "sol_prod = pd.concat([transp_cost, x_values_prod], axis=1)\n",
    "print(f\"This model has a total cost of {round(m_pwl2.ObjVal,2)}\")\n",
    "sol_prod[sol_prod.shipment > 0].drop('cost',axis=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Handling Uncertainty"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Simulation\n",
    "\n",
    "The first practical method discussed was basic Monte Carlo simulation. In the code above, we used [solution pools](https://docs.gurobi.com/projects/optimizer/en/current/features/solutionpool.html) to get the top two solutions. We already extracted the optimal solution. The code below queries the second best solution we asked Gurobi to find and puts the solutions together for comparison. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "This model has a total cost of 1755.44. The original cost was 1748.42\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th>cost</th>\n",
       "      <th>shipment</th>\n",
       "      <th>shipment2</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>production</th>\n",
       "      <th>distribution</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th rowspan=\"2\" valign=\"top\">Baltimore</th>\n",
       "      <th>Nashville</th>\n",
       "      <td>5.96</td>\n",
       "      <td>30.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Richmond</th>\n",
       "      <td>1.96</td>\n",
       "      <td>116.0</td>\n",
       "      <td>135.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">Cleveland</th>\n",
       "      <th>Columbia</th>\n",
       "      <td>2.43</td>\n",
       "      <td>89.0</td>\n",
       "      <td>89.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Indianapolis</th>\n",
       "      <td>2.37</td>\n",
       "      <td>95.0</td>\n",
       "      <td>40.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Nashville</th>\n",
       "      <td>4.13</td>\n",
       "      <td>0.0</td>\n",
       "      <td>30.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>St. Louis</th>\n",
       "      <td>4.88</td>\n",
       "      <td>0.0</td>\n",
       "      <td>41.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Little Rock</th>\n",
       "      <th>St. Louis</th>\n",
       "      <td>2.92</td>\n",
       "      <td>140.0</td>\n",
       "      <td>140.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Birmingham</th>\n",
       "      <th>Nashville</th>\n",
       "      <td>1.53</td>\n",
       "      <td>71.0</td>\n",
       "      <td>71.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"3\" valign=\"top\">Charleston</th>\n",
       "      <th>Indianapolis</th>\n",
       "      <td>2.61</td>\n",
       "      <td>0.0</td>\n",
       "      <td>55.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Lexington</th>\n",
       "      <td>1.61</td>\n",
       "      <td>121.0</td>\n",
       "      <td>121.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>St. Louis</th>\n",
       "      <td>4.60</td>\n",
       "      <td>41.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                          cost  shipment  shipment2\n",
       "production  distribution                           \n",
       "Baltimore   Nashville     5.96      30.0        0.0\n",
       "            Richmond      1.96     116.0      135.0\n",
       "Cleveland   Columbia      2.43      89.0       89.0\n",
       "            Indianapolis  2.37      95.0       40.0\n",
       "            Nashville     4.13       0.0       30.0\n",
       "            St. Louis     4.88       0.0       41.0\n",
       "Little Rock St. Louis     2.92     140.0      140.0\n",
       "Birmingham  Nashville     1.53      71.0       71.0\n",
       "Charleston  Indianapolis  2.61       0.0       55.0\n",
       "            Lexington     1.61     121.0      121.0\n",
       "            St. Louis     4.60      41.0        0.0"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m.Params.SolutionNumber = 1\n",
    "m_ObjVal2 = m.PoolObjVal\n",
    "print(f\"This model has a total cost of {round(m_ObjVal2,2)}. The original cost was {round(m_ObjVal,2)}\")\n",
    "x_values2 = pd.Series(m.getAttr('Xn', x), name = \"shipment2\", index = transp_cost.index)\n",
    "sol2 = pd.concat([transp_cost, x_values2], axis=1)\n",
    "sol_sim_comp = pd.concat([sol, sol2.drop('cost', axis = 1)], axis=1)\n",
    "sol_sim_comp[(sol_sim_comp.shipment > 0) | (sol_sim_comp.shipment2 > 0)]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fill in the blanks in the code below to simulate an uncertain transportation cost from Cleveland to Indianapolis. We will assume the cost is normally distributed with mean 2.3 and standard deviation 0.2. Generate 1000 samples. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0oAAAI1CAYAAADo2QpNAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy81sbWrAAAACXBIWXMAAA9hAAAPYQGoP6dpAAB4jUlEQVR4nO3deZxO9f//8eeF2XdmZYaxzIyxZpAs2ffI8qn4IKFNZIkolUwqisinlHZalFZSIrtkKUuyDSGMGMYwmDGLMXN+f/i5vs1lxjKuy5nlcb/drtt0neV1Xtc5RtfTOed9LIZhGAIAAAAAWJUyuwEAAAAAKGwISgAAAABgg6AEAAAAADYISgAAAABgg6AEAAAAADYISgAAAABgg6AEAAAAADYISgAAAABgg6AEAAAAADYISgBQyG3cuFH33nuvQkJC5OzsrODgYN1zzz3asGGDqX3NmTNHFotFhw4dcui6LVu2VMuWLa+5XHh4uLp06ZLnvM2bN8tisWjOnDkF6uHfJk2apAULFtzQOgCAooegBACF2JtvvqmmTZvqn3/+0ZQpU7R8+XK99tprOnr0qJo1a6aZM2ea3WKB3HXXXdqwYYNCQkKKXA8EJQAoGcqY3QAAIG/r1q3TyJEj1blzZ82fP19lyvzfX9m9e/dWjx49NGLECNWrV09Nmza9ZX2lp6fL1dX1pmoEBAQoICDATh0V3R4KIi0tTe7u7ma3AQDFHmeUAKCQmjx5siwWi2bNmpUrJElSmTJl9Pbbb8tiseiVV16RJC1YsEAWi0UrVqy4otasWbNksVi0fft267TNmzfr7rvvVtmyZeXq6qp69erpq6++yrXe5cvTli5dqkGDBikgIEDu7u7KzMzMs+dly5apW7duCg0Nlaurq6pVq6ZHH31USUlJedb992VvhmFoypQpqlSpklxdXRUTE6PFixff0D67EXn18Mcff6hLly4KDAyUi4uLypcvr7vuukv//POPJMlisej8+fP6+OOPZbFYZLFYcl0WuHPnTnXr1k1+fn5ydXXVbbfdpo8//viKbe/atUvt27eXu7u7AgICNHToUC1atEgWi0WrV6+2LteyZUvVqlVLv/zyi5o0aSJ3d3cNGjRIkvTll1+qffv2CgkJkZubm6Kjo/X000/r/PnzubY1YMAAeXp6as+ePerQoYM8PDwUEhJi/XOzceNGNWvWTB4eHoqMjLyi37S0ND355JOqXLmyXF1dVbZsWTVo0EBffPHFzex+ACj0OKMEAIVQdna2Vq1apQYNGig0NDTPZcLCwlS/fn2tXLlS2dnZ1i/4s2fPVps2bXItO2fOHMXExKhOnTqSpFWrVqljx45q1KiR3nnnHfn4+GjevHnq1auX0tLSNGDAgFzrDxo0SHfddZc+/fRTnT9/Xk5OTnn2dODAATVu3FgPPfSQfHx8dOjQIU2fPl3NmjXTjh078l1Pkl544QW98MILevDBB3XPPffoyJEjevjhh5Wdna2oqKjr2m+GYejixYtXTM/Ozr7muufPn1e7du1UuXJlvfXWWwoKCtLx48e1atUqpaSkSJI2bNig1q1bq1WrVho/frwkydvbW5K0d+9eNWnSRIGBgXrjjTdUrlw5ffbZZxowYIBOnDihsWPHSpISEhLUokULeXh4aNasWQoMDNQXX3yhxx9/PM++EhIS1K9fP40dO1aTJk1SqVKX/o1z37596ty5s0aOHCkPDw/t2bNHr776qn7//XetXLkyV42srCz17NlTgwcP1pgxY/T5559r3LhxOnfunL799ls99dRTCg0N1ZtvvqkBAwaoVq1aql+/viRp1KhR+vTTT/XSSy+pXr16On/+vHbu3KlTp05dzyEBgKLLAAAUOsePHzckGb17977qcr169TIkGSdOnDAMwzBGjRpluLm5GWfOnLEus3v3bkOS8eabb1qnVa9e3ahXr56RlZWVq16XLl2MkJAQIzs72zAMw5g9e7Yhyejfv/8V27487+DBg3n2lpOTY2RlZRmHDx82JBnff/99vusmJycbrq6uRo8ePXLVWLdunSHJaNGixVX3g2EYRqVKlQxJV33Nnj073x42b95sSDIWLFhw1e14eHgYDzzwwBXTe/fubbi4uBjx8fG5pnfq1Mlwd3e3HpMxY8YYFovF2LVrV67lOnToYEgyVq1aZZ3WokULQ5KxYsWKq/Z0eV+vWbPGkGT8+eef1nkPPPCAIcn49ttvrdOysrKMgIAAQ5KxdetW6/RTp04ZpUuXNkaNGmWdVqtWLaN79+5X3T4AFEdcegcARZhhGJIuXRImXTrzk56eri+//NK6zOzZs+Xi4qI+ffpIkvbv3689e/aob9++kqSLFy9aX507d1ZCQoL27t2bazv/+c9/rqufxMREDR48WGFhYSpTpoycnJxUqVIlSVJcXFy+623YsEEZGRnWni5r0qSJdf3r0axZM23atOmK1yeffHLNdatVqyY/Pz899dRTeuedd7R79+7r3q4krVy5Um3atFFYWFiu6QMGDFBaWpp1lMI1a9aoVq1aqlGjRq7l/vvf/+ZZ18/PT61bt75i+t9//60+ffooODhYpUuXlpOTk1q0aCHpyn1tsVjUuXNn6/syZcqoWrVqCgkJUb169azTy5Ytq8DAQB0+fNg67fbbb9fixYv19NNPa/Xq1UpPT7+e3QEARR6X3gFAIeTv7y93d3cdPHjwqssdOnRI7u7uKlu2rCSpZs2aatiwoWbPnq1HHnlE2dnZ+uyzz9StWzfrMidOnJAkPfnkk3ryySfzrGt7T9H1jAyXk5Oj9u3b69ixYxo/frxq164tDw8P5eTk6I477rjqF+zLl3EFBwdfMS+vafnx8fFRgwYNrnt523XXrFmjl19+Wc8884ySk5MVEhKihx9+WM8999xVLxuULn2GvPZT+fLlrfMv/6xcufIVywUFBeVZN6+aqampuvPOO+Xq6qqXXnpJkZGRcnd315EjR9SzZ88r9rW7u/sVA3A4Oztb/0zYTs/IyLC+f+ONNxQaGqovv/xSr776qlxdXdWhQwdNnTpVERERefYMAMUBQQkACqHSpUurVatWWrJkif75558871P6559/tGXLFnXq1EmlS5e2Th84cKCGDBmiuLg4/f3330pISNDAgQOt8/39/SVJ48aNU8+ePfPcvu09QZfPWF3Nzp079eeff2rOnDl64IEHrNP3799/zXXLlSsnSTp+/PgV844fP67w8PBr1rCH2rVra968eTIMQ9u3b9ecOXM0ceJEubm56emnn77quuXKlVNCQsIV048dOybp//Z7uXLlrGH13/L67FLe+37lypU6duyYVq9ebT2LJElnzpy5ao8F4eHhYb1/7MSJE9azS127dtWePXvsvj0AKCy49A4ACqlx48bJMAwNGTLkisEIsrOz9dhjj8kwDI0bNy7XvP/+979ydXXVnDlzNGfOHFWoUEHt27e3zo+KilJERIT+/PNPNWjQIM+Xl5fXDfd7+Qu9i4tLrunvvvvuNde944475Orqqrlz5+aavn79+lyXgd0qFotFdevW1euvvy5fX19t3brVOs/FxSXPs2Nt2rSxBph/++STT+Tu7q477rhDktSiRQvt3Lnzikv75s2bd0P9Xe7l365nX9+MoKAgDRgwQP/973+1d+9epaWlOXR7AGAmzigBQCHVtGlTzZgxQyNHjlSzZs30+OOPq2LFioqPj9dbb72l3377TTNmzFCTJk1yrefr66sePXpozpw5OnPmjJ588knrSGmXvfvuu+rUqZM6dOigAQMGqEKFCjp9+rTi4uK0detWff311zfcb/Xq1VW1alU9/fTTMgxDZcuW1Q8//KBly5Zdc10/Pz89+eSTeumll/TQQw/p3nvv1ZEjRxQbG3tDl97djB9//FFvv/22unfvripVqsgwDH333Xc6c+aM2rVrZ12udu3aWr16tX744QeFhITIy8tLUVFRmjBhgn788Ue1atVKzz//vMqWLau5c+dq0aJFmjJlinx8fCRJI0eO1EcffaROnTpp4sSJCgoK0ueff249O2N7rPLSpEkT+fn5afDgwZowYYKcnJw0d+5c/fnnn3bfL40aNVKXLl1Up04d+fn5KS4uTp9++qkaN27M85wAFGucUQKAQmzYsGFat26dQkNDNXr0aLVu3VqjRo1SSEiIfv31Vw0bNizP9QYOHKjExERduHDhiqG+JalVq1b6/fff5evrq5EjR6pt27Z67LHHtHz5crVt27ZAvTo5OemHH35QZGSkHn30Uf33v/9VYmKili9ffl3rT5w4UZMnT9bSpUt19913680339Q777xz3UOD36yIiAj5+vpqypQpuvvuu3Xvvfdq69atmjNnjh5++GHrcv/73/8UERGh3r17q2HDhnr00UclXTpTt379ekVFRWno0KHq3r27du7cqdmzZ2vMmDHW9cuXL681a9YoMjJSgwcPVt++feXs7KyJEydKuhR0r6VcuXJatGiR3N3d1a9fPw0aNEienp65BvGwl9atW2vhwoUaOHCg2rdvrylTpqh///764Ycf7L4tAChMLMblIZMAAIBpHnnkEX3xxRc6deqUnJ2dzW4HAEo8Lr0DAOAWmzhxosqXL68qVaooNTVVP/74oz744AM999xzhCQAKCQISgAA3GJOTk6aOnWq/vnnH128eFERERGaPn26RowYYXZrAID/j0vvAAAAAMAGgzkAAAAAgA2CEgAAAADYICgBAAAAgI1iP5hDTk6Ojh07Ji8vL+uTzAEAAACUPIZhKCUlReXLl7/mA76LfVA6duyYwsLCzG4DAAAAQCFx5MgRhYaGXnWZYh+UvLy8JF3aGd7e3iZ3AwAlU/XqUkKCFBIi7dljdjcm+6G6lJ4guYVIXUv6zgCAW+vcuXMKCwuzZoSrKfZB6fLldt7e3gQlADDJ5asbSpWSSvxfxe6lJIskN3YGAJjlem7JYTAHAAAAALBBUAIAAAAAGwQlAAAAALBR7O9RAgCYb9MmKTtbKl3a7E4KgY6bJCNbsrAzYH/Z2dnKysoyuw3ANE5OTiptp//ZEJQAAA4XEmJ2B4WIGzsD9mcYho4fP64zZ86Y3QpgOl9fXwUHB9/0M1QJSgAAAEXc5ZAUGBgod3f3m/6CCBRFhmEoLS1NiYmJkqSQm/xXOoISAABAEZadnW0NSeXKlTO7HcBUbm5ukqTExEQFBgbe1GV4BCUAgMO9956Umip5ekqPPGJ2Nybb/56UlSo5eUrVSvrOgD1cvifJ3d3d5E6AwuHy70JWVhZBCQBQuE2cKB09KlWoQFDSjolS+lHJrQJBCXbF5XbAJfb6XWB4cAAAAACwwRklAACAYio+Pl5JSUm3ZFv+/v6qWLHiLdnWjYiNjdWCBQu0bdu2YrEdi8Wi+fPnq3v37jdVJzw8XCNHjtTIkSPt0ldxRFACAAAohuLj41W9erTS09Nuyfbc3Ny1Z0/cDYelI0eOKDY2VosXL1ZSUpJCQkLUvXt3Pf/88zc8OEVeIeLJJ5/UsGHDbqiOo3z77beaMmWK9uzZo5ycHFWsWFEdO3bUtGnTHLbNOXPmaOTIkVcMHb9p0yZ5eHg4bLvFAUEJAACgGEpKSlJ6epp69PhMAQHRDt3WyZNxmj+/n5KSkm4oKP39999q3LixIiMj9cUXX6hy5cratWuXxowZo8WLF2vjxo0qW7bsTfXm6ekpT0/Pm6phD8uXL1fv3r01adIk3X333bJYLNq9e7dWrFhhSj8BAQGmbLco4R4lAACAYiwgIFohITEOfRU0iA0dOlTOzs5aunSpWrRooYoVK6pTp05avny5jh49qmeffda6bHh4uF588UX16dNHnp6eKl++vN58881c8yWpR48eslgs1vexsbG67bbbrMsNGDBA3bt316RJkxQUFCRfX1+98MILunjxosaMGaOyZcsqNDRUH330Ua5en3rqKUVGRsrd3V1VqlTR+PHjrSMOXo8ff/xRzZo105gxYxQVFaXIyEh1794912eQpFmzZqlq1apydnZWVFSUPv3003xrrl69WhaLJdfZom3btslisejQoUNavXq1Bg4cqLNnz8pischisSg2Nta6v2bMmGFdLz4+Xt26dZOnp6e8vb1133336cSJE9b5l/fjp59+qvDwcPn4+Kh3795KSUm57n1Q1BCUAAAAcMudPn1aP//8s4YMGWJ99s1lwcHB6tu3r7788ksZhmGdPnXqVNWpU0dbt27VuHHj9MQTT2jZsmWSLl1KJkmzZ89WQkKC9X1eVq5cqWPHjumXX37R9OnTFRsbqy5dusjPz0+//fabBg8erMGDB+vIkSPWdby8vDRnzhzt3r1b//vf//T+++/r9ddfv+7PGxwcrF27dmnnzp35LjN//nyNGDFCo0eP1s6dO/Xoo49q4MCBWrVq1XVv59+aNGmiGTNmyNvbWwkJCUpISNCTTz55xXKGYah79+46ffq01qxZo2XLlunAgQPq1atXruUOHDigBQsW6Mcff9SPP/6oNWvW6JVXXilQb0UBQQkAAAC33L59+2QYhqKj8z4bFR0dreTkZJ08edI6rWnTpnr66acVGRmpYcOG6Z577rGGlcuXkvn6+io4OPiql5aVLVtWb7zxhqKiojRo0CBFRUUpLS1NzzzzjCIiIjRu3Dg5Oztr3bp11nWee+45NWnSROHh4eratatGjx6tr7766ro/77Bhw9SwYUPVrl1b4eHh6t27tz766CNlZmZal3nttdc0YMAADRkyRJGRkRo1apR69uyp11577bq382/Ozs7y8fGRxWJRcHCwgoOD87wMcfny5dq+fbs+//xz1a9fX40aNdKnn36qNWvW5AqcOTk5mjNnjmrVqqU777xT999/v2mXDt4KBCUAAAAUOpfPJP37mTiNGzfOtUzjxo0VFxd3w7Vr1qypUqX+72twUFCQateubX1funRplStXTomJidZp33zzjZo1a2YNG+PHj1d8fPx1b9PDw0OLFi3S/v379dxzz8nT01OjR4/W7bffrrS0SwNuxMXFqWnTprnWa9q0aYE+442Ii4tTWFiYwsLCrNNq1KghX1/fXNsODw+Xl5eX9X1ISEiufVTcEJQAAA4XGSnVqHHpZ4nnHSn51Lj0EyjBqlWrZh3QIC979uyRn5+f/P39r1qnIA8XdXJyuqJGXtNycnIkSRs3blTv3r3VqVMn/fjjj/rjjz/07LPP6sKFCze87apVq+qhhx7SBx98oK1bt2r37t368ssv8/08hmHk+xkvh71/X554I/dNXWsbttOvto+KI0a9AwA43MqVZndw4xz2/Bm/1wrt82aAW6lcuXJq166d3n77bT3xxBO57lM6fvy45s6dq/79++f6or5x48ZcNTZu3Kjq1atb3zs5OSk7O9vuva5bt06VKlXKNbjE4cOHb7pueHi43N3ddf78eUmXLjf89ddf1b9/f+sy69evz/fyxMuXFyYkJMjPz0+SrniOk7Oz8zX3SY0aNRQfH68jR45Yzyrt3r1bZ8+ezXfbJQFBCQAAG/Hx8YqOrq60tHSH1Hd3d1Nc3B7CEkq8mTNnqkmTJurQoYNeeumlXMODV6hQQS+//HKu5detW6cpU6aoe/fuWrZsmb7++mstWrTIOj88PFwrVqxQ06ZN5eLiYg0PN6tatWqKj4/XvHnz1LBhQy1atEjz58+/oRqxsbFKS0tT586dValSJZ05c0ZvvPGGsrKy1K5dO0nSmDFjdN999ykmJkZt2rTRDz/8oO+++07Lly/Pt6+wsDDFxsbqpZde0r59+654JlN4eLhSU1O1YsUK1a1bV+7u7nJ3d8+1TNu2bVWnTh317dtXM2bM0MWLFzVkyBC1aNFCDRo0uKHPWZwQlAAAsJGUlKS0tHR9NqmHoqvY91kjcX+fVL9n5t/w82aAgjp50rH3t9zMNiIiIrR582bFxsaqV69eOnXqlIKDg9W9e3dNmDDhimcojR49Wlu2bNELL7wgLy8vTZs2TR06dLDOnzZtmkaNGqX3339fFSpU0KFDh27mY1l169ZNTzzxhB5//HFlZmbqrrvu0vjx461DbV+PFi1a6K233lL//v114sQJ+fn5qV69elq6dKmioqIkSd27d9f//vc/TZ06VcOHD1flypU1e/ZstWzZMs+aTk5O+uKLL/TYY4+pbt26atiwoV566SXde++91mWaNGmiwYMHW/fvhAkTrujbYrFowYIFGjZsmJo3b65SpUqpY8eOVwxdXtJYjH9f1FgMnTt3Tj4+Pjp79qy8vb3NbgcAUARs3bpV9evX15Z5jygmOsS+teMSVL/3e9qyZYtiYmLsWhslU0ZGhg4ePKjKlSvL1dXVOj0+Pl7Vq0crPT3tlvTh5uauPXviHPYPAOHh4Ro5cqRGjhzpkPooPvL7nZBuLBtwRgkA4HB9+0pJSZK/vzR3rtndmCs8Z6WWPCWFH3tWillsdjsoxipWrKg9e+Icc69dHrj3DsUNQQkA4HBr1khHj0oVKpjdifk8jQR1qCNdSNtqdisoASpWrEh4AQqIoAQAAIBCz173GwHXi+coAQAAAIANghIAAAAA2CAoAQAAAIANghIAAAAA2CAoAQAAAIANghIAAAAA2GB4cAAAgGIqPj6eB846yIABA3TmzBktWLDA7FZuGYvFovnz56t79+43VSc8PFwjR47UyJEj7dKXoxCUAAAO9/DD0tmzko+P2Z2YL8lSXfMWbVHv/j1U3uxmUKzFx8crOrq60tLSb8n23N3dFBe354bCUmJiosaPH6/FixfrxIkT8vPzU926dRUbG6vGjRs7sFvHW716tVq1amV97+rqqipVqmjEiBF65JFH7LKNQ4cOqXLlyvrjjz902223XXXZb7/9VlOmTNGePXuUk5OjihUrqmPHjpo2bZpdesnLnDlzNHLkSJ05cybX9E2bNsnDw8Nh27UXghIAwOEmTDC7g8LjeKn6Gj13i1qOeoSgBIdKSkpSWlq6PpvUQ9FVAhy6rbi/T6rfM/OVlJR0Q0HpP//5j7KysvTxxx+rSpUqOnHihFasWKHTp087sNtba+/evfL29lZ6erp++OEHPfbYY6pataratGlzy3pYvny5evfurUmTJunuu++WxWLR7t27tWLFilvWw78FBDj2z6O9cI8SAABAMRZdJUAx0SEOfRUkiJ05c0a//vqrXn31VbVq1UqVKlXS7bffrnHjxumuu+6yLnf27Fk98sgjCgwMlLe3t1q3bq0///wzV62FCxeqQYMGcnV1lb+/v3r27Gmdl5ycrP79+8vPz0/u7u7q1KmT9u3bZ50/Z84c+fr66ueff1Z0dLQ8PT3VsWNHJSQkWJfJzs7WqFGj5Ovrq3Llymns2LEyDOO6PmdgYKCCg4NVuXJlDR8+XOHh4dq6dat1vmEYmjJliqpUqSI3NzfVrVtX33zzTa7++/btq4CAALm5uSkiIkKzZ8+WJFWuXFmSVK9ePVksFrVs2TLPHn788Uc1a9ZMY8aMUVRUlCIjI9W9e3e9+eabuZabNWuWqlatKmdnZ0VFRenTTz/N93OtXr1aFosl19mibdu2yWKx6NChQ1q9erUGDhyos2fPymKxyGKxKDY2VtKlS+9mzJhhXS8+Pl7dunWTp6envL29dd999+nEiRPW+bGxsbrtttv06aefKjw8XD4+Purdu7dSUlKuuu9vFkEJAAAAt5ynp6c8PT21YMECZWZm5rmMYRi66667dPz4cf3000/asmWLYmJi1KZNG+tZp0WLFqlnz56666679Mcff2jFihVq0KCBtcaAAQO0efNmLVy4UBs2bJBhGOrcubOysrKsy6Slpem1117Tp59+ql9++UXx8fF68sknrfOnTZumjz76SB9++KF+/fVXnT59WvPnz7+hz2sYhpYsWaIjR46oUaNG1unPPfecZs+erVmzZmnXrl164okn1K9fP61Zs0aSNH78eO3evVuLFy9WXFycZs2aJX9/f0nS77//LunSGaOEhAR99913eW47ODhYu3bt0s6dO/Ptb/78+RoxYoRGjx6tnTt36tFHH9XAgQO1atWqG/qclzVp0kQzZsyQt7e3EhISlJCQkGuf/nu/dO/eXadPn9aaNWu0bNkyHThwQL169cq13IEDB7RgwQL9+OOP+vHHH7VmzRq98sorBertenHpHQAAAG65MmXKaM6cOXr44Yf1zjvvKCYmRi1atFDv3r1Vp04dSdKqVau0Y8cOJSYmysXFRZL02muvacGCBfrmm2/0yCOP6OWXX1bv3r31wgsvWGvXrVtXkrRv3z4tXLhQ69atU5MmTSRJc+fOVVhYmBYsWKB7771XkpSVlaV33nlHVatWlSQ9/vjjmjhxorXejBkzNG7cOP3nP/+RJL3zzjv6+eefr+tzhoaGSpIyMzOVk5OjiRMnqnnz5pKk8+fPa/r06Vq5cqX1nqwqVaro119/1bvvvqsWLVooPj5e9erVs4a/8PBwa+3Ll7CVK1dOwcHB+fYwbNgwrV27VrVr11alSpV0xx13qH379urbt2+u/TpgwAANGTJEkjRq1Cht3LhRr732Wq57ra6Xs7OzfHx8ZLFYrtrb8uXLtX37dh08eFBhYWGSpE8//VQ1a9bUpk2b1LBhQ0lSTk6O5syZIy8vL0nS/fffrxUrVujll1++4d6uF2eUAAAOFxoqWSyXfpZ0tbLnypgr1drfyexWANP95z//0bFjx7Rw4UJ16NBBq1evVkxMjObMmSNJ2rJli1JTU1WuXDnrGShPT08dPHhQBw4ckHTpcq/87veJi4tTmTJlcp3BKVeunKKiohQXF2ed5u7ubg1JkhQSEqLExERJly79S0hIyDW4RJkyZXKdtbqatWvXatu2bdq2bZs++OADTZo0SbNmzZIk7d69WxkZGWrXrl2uz/fJJ59YP99jjz2mefPm6bbbbtPYsWO1fv3669ruv3l4eGjRokXav3+/nnvuOXl6emr06NG6/fbblZaWZt1XTZs2zbVe06ZNc+0nR4iLi1NYWJg1JElSjRo15Ovrm2vb4eHh1pAk5T5GjsIZJQAAAJjG1dVV7dq1U7t27fT888/roYce0oQJEzRgwADl5OQoJCREq1evvmI9X19fSZKbm1u+tfO7j8gwDFksFut7JyenXPMtFst134N0LZUrV7b2WrNmTf322296+eWX9dhjjyknJ0fSpcsHK1SokGu9y2d6OnXqpMOHD2vRokVavny52rRpo6FDh+q111674V6qVq2qqlWr6qGHHtKzzz6ryMhIffnllxo4cKAk5don0pX76d9KlSplXeayf1/OeL3y28b1HKPL+89ROKMEAACAQqNGjRo6f/68JCkmJkbHjx9XmTJlVK1atVyvy/fp1KlTJ9/R22rUqKGLFy/qt99+s047deqU/vrrL0VHR19XPz4+PgoJCdHGjRut0y5evKgtW7YU6POVLl1a6enp1v5cXFwUHx9/xef79xmWgIAADRgwQJ999plmzJih9957T9Kly9ukS4NN3Kjw8HC5u7tb93V0dLR+/fXXXMusX78+3/10+bK/fw96sW3btlzLODs7X7O3GjVqKD4+XkeOHLFO2717t86ePXvdx8hROKMEAABQjMX9fbJQbuPUqVO69957NWjQINWpU0deXl7avHmzpkyZom7dukmS2rZtq8aNG6t79+569dVXFRUVpWPHjumnn35S9+7d1aBBA02YMEFt2rRR1apV1bt3b128eFGLFy/W2LFjFRERoW7duunhhx/Wu+++Ky8vLz399NOqUKGCdRvXY8SIEXrllVcUERGh6OhoTZ8+/YpnA+UnMTFRGRkZyszM1O+//65PP/1U99xzjyTJy8tLTz75pJ544gnl5OSoWbNmOnfunNavXy9PT0898MADev7551W/fn3VrFlTmZmZ+vHHH60BIjAwUG5ublqyZIlCQ0Pl6uoqnzweWBcbG6u0tDR17txZlSpV0pkzZ/TGG28oKytL7dq1kySNGTNG9913n3WwjB9++EHfffedli9fnufnuhzmYmNj9dJLL2nfvn1XPJMpPDxcqampWrFiherWrSt3d3e5u7vnWqZt27aqU6eO+vbtqxkzZujixYsaMmSIWrRocd2XNzoKQQkAAKAY8vf3l7u7m/o9c2OjsxWUu7ub9SzP9fD09FSjRo30+uuv68CBA8rKylJYWJgefvhhPfPMM5IuXV71008/6dlnn9WgQYN08uRJBQcHq3nz5goKCpIktWzZUl9//bVefPFFvfLKK/L29rYOliBJs2fP1ogRI9SlSxdduHBBzZs3108//XTFpVxXM3r0aCUkJGjAgAEqVaqUBg0apB49eujs2bPXXDcqKkrSpfuawsLC9Oijj1qHyZakF198UYGBgZo8ebL+/vtv+fr6KiYmxroPnJ2dNW7cOB06dEhubm668847NW/ePGvNN954QxMnTtTzzz+vO++8M8/LFFu0aKG33npL/fv3tz7Yt169elq6dKm1v+7du+t///ufpk6dquHDh6ty5cqaPXt2vkOOOzk56YsvvtBjjz2munXrqmHDhnrppZesA2RIl0a+Gzx4sHr16qVTp05pwoQJuT67dOkYL1iwQMOGDVPz5s1VqlQpdezY8Yqhy81gMex1AWYhde7cOfn4+Ojs2bPy9vY2ux0AKJFCQ6WjR6UKFaR//jG7m2vbunWr6tevry3zHlFMdIhda1/Y+ZqcdV4XygTK+b4T114BuIaMjAwdPHhQlStXlqura6558fHxSkpKuiV9+Pv739DDZgFHudrvxI1kA84oAQAAFFMVK1YkvAAFxGAOAAAAAGCDoAQAAAAANghKAAAAAGCDoAQAcLjPPpOWLLn0s6Q7VKqVOrwiHQp50exWUMwU8/G5gOtmr98FBnMAADhcPqPLlkiplvJaukNK9TD3+SAoPi4Pc52WliY3NzeTuwHMl5aWJkk3NAR8XghKAAAARVjp0qXl6+urxMRESZK7u7ssFovJXQG3nmEYSktLU2Jionx9fVW6dOmbqkdQAgAAKOKCg4MlyRqWgJLM19fX+jtxMwhKAACHW71aysyUXFy4DM/TOKb2tSXP85slxZjdDooJi8WikJAQBQYGKisry+x2ANM4OTnd9JmkywhKAACH69dPOnpUqlBB+ucfs7sxV3jOKv38tHQhYbykR8xuB8VM6dKl7fYlESjpGPUOAAAAAGxwRgkAABNczMrSzq1b7VrT399fFStWtGtNACipCEoAANxC2Tk5UinpdHKy6tevb9fa7u5uiovbQ1gCADsgKAEAcAsZOYZUSvL1dNGWeQ/YrW7c3yfV75n5SkpKIigBgB0QlAAAMEGZMhbFRIeY3QYAIB8M5gAAAAAANghKAAAAAGCDoAQAAAAANghKAAAAAGCDoAQAcLh//pEM49LPkm7+sc6y9JXWnOlsdisAgKsgKAEAAACADYISAAAAANggKAEAAACADR44CwBwuBdekM6elXx8pAkTzO7GXLW94zStr1TVNU5SbbPbAQDkg6AEAHC499+Xjh6VKlQgKFXzPKg6naWMnINmtwIAuAouvQMAAAAAGwQlAAAAALBBUAIAAAAAG9yjBAAosuLj45WUlGT3unFxcXavCQAoWghKAIAiKT4+XtHR1ZWWlu6wbaSkpDqsNgCgcCMoAQCKpKSkJKWlpeuzST0UXSXArrV/WrtP499apYyMDLvWBQAUHQQlAECRFl0lQDHRIXatGXfQ/pfzAQCKFgZzAAAAAAAbnFECADhcixZSUpLk7292J+Y7kRmgv3YfUe2oAJU3uxkAQL4ISgAAh5s71+wOCo/1pxqq36tHtOTNhgQlACjEuPQOAAAAAGwUmqA0efJkWSwWjRw50jrNMAzFxsaqfPnycnNzU8uWLbVr1y7zmgQAAABQIhSKoLRp0ya99957qlOnTq7pU6ZM0fTp0zVz5kxt2rRJwcHBateunVJSUkzqFAAAAEBJYHpQSk1NVd++ffX+++/Lz8/POt0wDM2YMUPPPvusevbsqVq1aunjjz9WWlqaPv/8cxM7BgDcqNatpZo1L/0s6doErtXOV6UGXmvNbgUAcBWmB6WhQ4fqrrvuUtu2bXNNP3jwoI4fP6727dtbp7m4uKhFixZav379rW4TAHAT/vpL2r370s+SzrtMimqGSh6luToCAAozU0e9mzdvnrZu3apNmzZdMe/48eOSpKCgoFzTg4KCdPjw4XxrZmZmKjMz0/r+3LlzduoWAAAAQElh2hmlI0eOaMSIEfrss8/k6uqa73IWiyXXe8Mwrpj2b5MnT5aPj4/1FRYWZreeAQAAAJQMpgWlLVu2KDExUfXr11eZMmVUpkwZrVmzRm+88YbKlCljPZN0+czSZYmJiVecZfq3cePG6ezZs9bXkSNHHPo5AAAAABQ/pl1616ZNG+3YsSPXtIEDB6p69ep66qmnVKVKFQUHB2vZsmWqV6+eJOnChQtas2aNXn311Xzruri4yMXFxaG9AwAAACjeTAtKXl5eqlWrVq5pHh4eKleunHX6yJEjNWnSJEVERCgiIkKTJk2Su7u7+vTpY0bLAAAAAEoIUwdzuJaxY8cqPT1dQ4YMUXJysho1aqSlS5fKy8vL7NYAAAAAFGOFKiitXr0613uLxaLY2FjFxsaa0g8AAACAksn05ygBAAAAQGFTqM4oAQCKp+efl1JTJU9Pszsx346z0fpy8R96sEe0aprdDAAgXwQlAIDDPfKI2R0UHvvPV9bri/9Qh86VCUoAUIhx6R0AAAAA2CAoAQAAAIANLr0DADhcQoKUnS2VLi2FhJjdjblcS6WrQlnJ2ZJudisAgKvgjBIAwOEaNpTCwi79LOk6Ba/SP29KjX1Wmd0KAOAqCEoAAAAAYIOgBAAAAAA2CEoAAAAAYIOgBAAAAAA2CEoAAAAAYIOgBAAAAAA2CEoAAAAAYIOgBAAAAAA2CEoAAAAAYIOgBABwuBUrpJ07L/0s6ZYn3qmaY6VN5+40uxUAwFWUMbsBAEDxFxVldgeFR8pFL+0+KqXleJndCgDgKjijBAAAAAA2CEoAAAAAYINL7wAADvf551JamuTuLvXpY3Y35gp3P6IHW0ohzkck1Ta7HQBAPghKAACHGztWOnpUqlCBoFTPd4eaPixl5OyQ1NnsdgAA+eDSOwAAAACwQVACAAAAABsEJQAAAACwQVACAAAAABsEJQAAAACwQVACAAAAABsEJQAAAACwQVACAAAAABs8cBYA4HDBwbl/lmTp2a46fS5DXt6ucjW7GQBAvjijBABwuM2bpX/+ufSzpFtyorXChkkbz7U2uxUAwFUQlAAAAADABkEJAAAAAGwQlAAAAADABoM5AAAc7tFHpdOnpbJlpXffNbsbc93ut1VfDZdquG+VVNvsdgAA+SAoAQAcbtEi6ehRqUIFszsxXwW344poJGXkHDe7FQDAVXDpHQAAAADYICgBAAAAgA2CEgAAAADYICgBAAAAgA2CEgAAAADYICgBAAAAgA2CEgAAAADYICgBAAAAgA0eOAsAcLj//ldKTpb8/MzuxHyH0sK0fss+tbkjTJXNbgYAkC+CEgDA4aZONbuDwuOPM7X18Af7tKRubYISABRiXHoHAAAAADYISgAAAABgg6AEAAAAADYISgAAh6teXfL2vvSzpOsSslRnP5Ca+iw1uxUAwFUQlAAADpeaKqWkXPpZ0jlZLsrbTSpjuWh2KwCAqyAoAQAAAIANhgcHAKAYiYuLc0hdf39/VaxY0SG1AaAwIigBAFAMJCSlyiKpX79+Dqnv7u6muLg9hCUAJQZBCQCAYuBMSoYMSTPHtlLjmAi71o77+6T6PTNfSUlJBCUAJQZBCQCAYqRamJ9iokPMbgMAijwGcwAAAAAAGwQlAAAAALBBUAIAAAAAG9yjBABwuHfekdLTJTc3szsx3++n6+mtLzZo3KB6ijG7GQBAvghKAACH69LF7A4Kj6MZIfrmd+mh+xlwAQAKMy69AwAAAAAbnFECADhUfHy8kpKS7F43Li7O7jUBALiMoAQAcJj4+HhFR1dXWlq0JGdJFyRttes2UlJS7VrP0co6JeuOapJ36WSzWwEAXAVBCQDgMElJSUpLS5ef91IlnyunwLLJWvz203ap/dPafRr/1iplZGTYpd6t0iJggza8IGXkbJDU3Ox2AAD5ICgBABzOqUyp//+ztGKi7TOIQdxB+1/OBwDAZQzmAAAAAAA2CEoAAAAAYIOgBAAAAAA2CEoAAAAAYIOgBAAAAAA2CEoAAAAAYIOgBAAAAAA2CEoAAAAAYIOgBABwuG+mTdDZdZMVt2Cm2a2Y7oeEdvJ+SPr1TDuzWwEAXEUZsxsAABR/Hm6Z8vbMNLuNQuGi4aSUdClbTma3AgC4Cs4oAQAAAIANghIAAAAA2ODSOwCAw322qK0Wrg6Qt0emRvXfYHY7pqrutU8TekqVXPdJqm12OwCAfBCUAAAON3dRWyWe9lOFwHMlPihFe+1T/f9IGTn7zG4FAHAVXHoHAAAAADYISgAAAABgw9SgNGvWLNWpU0fe3t7y9vZW48aNtXjxYut8wzAUGxur8uXLy83NTS1bttSuXbtM7BgAAABASWBqUAoNDdUrr7yizZs3a/PmzWrdurW6detmDUNTpkzR9OnTNXPmTG3atEnBwcFq166dUlJSzGwbAAAAQDFnalDq2rWrOnfurMjISEVGRurll1+Wp6enNm7cKMMwNGPGDD377LPq2bOnatWqpY8//lhpaWn6/PPPzWwbAAAAQDFXaO5Rys7O1rx583T+/Hk1btxYBw8e1PHjx9W+fXvrMi4uLmrRooXWr1+fb53MzEydO3cu1wsAAAAAboTpQWnHjh3y9PSUi4uLBg8erPnz56tGjRo6fvy4JCkoKCjX8kFBQdZ5eZk8ebJ8fHysr7CwMIf2DwAAAKD4MT0oRUVFadu2bdq4caMee+wxPfDAA9q9e7d1vsViybW8YRhXTPu3cePG6ezZs9bXkSNHHNY7AAAAgOLJ9AfOOjs7q1q1apKkBg0aaNOmTfrf//6np556SpJ0/PhxhYSEWJdPTEy84izTv7m4uMjFxcWxTQMAbkj1yvGqEpqqAL80s1sx3ekLvvrz4HFVreQrV7ObAQDky/QzSrYMw1BmZqYqV66s4OBgLVu2zDrvwoULWrNmjZo0aWJihwCAG/X6mLe14dMPtfCNL8xuxXRrkpqoSaz0Ryr/LwOAwszUM0rPPPOMOnXqpLCwMKWkpGjevHlavXq1lixZIovFopEjR2rSpEmKiIhQRESEJk2aJHd3d/Xp08fMtgEAAAAUc6YGpRMnTuj+++9XQkKCfHx8VKdOHS1ZskTt2rWTJI0dO1bp6ekaMmSIkpOT1ahRIy1dulReXl5mtg0AAACgmDM1KH344YdXnW+xWBQbG6vY2Nhb0xAAAAAAqBAM5gAAKP6emDpEF7LKKsAvrcTfp9TCf73Wx0pVPddLqm12OwCAfBCUAAAOt+dgRSWe9lOFQB4CXtb5jEIjpIycM2a3AgC4ikI36h0AAAAAmI2gBAAAAAA2CEoAAAAAYIOgBAAAAAA2CEoAAAAAYIOgBAAAAAA2CEoAAAAAYIPnKAEAgOsSFxfnkLr+/v6qWLGiQ2oDQEERlAAADtf3ruXy9giQt0em2a2YLi4lQj+s3KHenSJU3exmrlNCUqoskvr16+eQ+u7uboqL20NYAlCoFCgoHTx4UJUrV7Z3LwCAYqrfXcsVEx1idhuFwp6UCL3w3Q41blV0gtKZlAwZkmaObaXGMRF2rR3390n1e2a+kpKSCEoACpUCBaVq1aqpefPmevDBB3XPPffI1dXV3n0BAIBCplqYH4EXQIlRoMEc/vzzT9WrV0+jR49WcHCwHn30Uf3+++/27g0AAAAATFGgoFSrVi1Nnz5dR48e1ezZs3X8+HE1a9ZMNWvW1PTp03Xy5El79wkAKMLOp7voXKqLUs47m92K6cpYsuTlJpVWltmtAACu4qaGBy9Tpox69Oihr776Sq+++qoOHDigJ598UqGhoerfv78SEhLs1ScAoAi7Z/QL8mk6TtHdHze7FdN1DVmmcx9IzXyXmd0KAOAqbioobd68WUOGDFFISIimT5+uJ598UgcOHNDKlSt19OhRdevWzV59AgAAAMAtU6DBHKZPn67Zs2dr79696ty5sz755BN17txZpUpdyl2VK1fWu+++q+rVi8p4PgAAAADwfwoUlGbNmqVBgwZp4MCBCg4OznOZihUr6sMPP7yp5gAAAADADAUKSvv27bvmMs7OznrggQcKUh4AAAAATFWge5Rmz56tr7/++orpX3/9tT7++OObbgoAAAAAzFSgoPTKK6/I39//iumBgYGaNGnSTTcFAAAAAGYqUFA6fPiwKleufMX0SpUqKT4+/qabAgAAAAAzFSgoBQYGavv27VdM//PPP1WuXLmbbgoAAAAAzFSgoNS7d28NHz5cq1atUnZ2trKzs7Vy5UqNGDFCvXv3tnePAAAAAHBLFWjUu5deekmHDx9WmzZtVKbMpRI5OTnq378/9ygBAK4w/cm3VSU0SM5O2Wa3Yro1Jxtr4jurNH1UYzU2uxkAQL4KFJScnZ315Zdf6sUXX9Sff/4pNzc31a5dW5UqVbJ3fwCAYiC6SrxiorPMbqNQOJ3lp437pXPZfma3AgC4igIFpcsiIyMVGRlpr14AAAAAoFAoUFDKzs7WnDlztGLFCiUmJionJyfX/JUrV9qlOQAAAAAwQ4GC0ogRIzRnzhzdddddqlWrliwWi737AgAUI79sqa0DRwLl5nJRXVr8ZXY7pqrgmqB7bpcCnBIk1Ta7HQBAPgoUlObNm6evvvpKnTt3tnc/AIBiaPKHfZV42k8VAs/pnxbTzW7HVLeX/UMtR0gZOX9Iam92OwCAfBRoeHBnZ2dVq1bN3r0AAAAAQKFQoKA0evRo/e9//5NhGPbuBwAAAABMV6BL73799VetWrVKixcvVs2aNeXk5JRr/nfffWeX5gAAAADADAUKSr6+vurRo4e9ewEAAACAQqFAQWn27Nn27gMAAAAACo0C3aMkSRcvXtTy5cv17rvvKiUlRZJ07Ngxpaam2q05AAAAADBDgc4oHT58WB07dlR8fLwyMzPVrl07eXl5acqUKcrIyNA777xj7z4BAAAA4JYp0BmlESNGqEGDBkpOTpabm5t1eo8ePbRixQq7NQcAAAAAZijwqHfr1q2Ts7NzrumVKlXS0aNH7dIYAKD4cHPNlJdHpjzdL5jdiumyjDI6ly6Vci7Q/4IBALdIgc4o5eTkKDs7+4rp//zzj7y8vG66KQBA8fLd9Ak6t36y9nw/0+xWTPdjQnv5PCStO9ve7FYAAFdRoKDUrl07zZgxw/reYrEoNTVVEyZMUOfOne3VGwAAAACYokDn/V9//XW1atVKNWrUUEZGhvr06aN9+/bJ399fX3zxhb17BAAAAIBbqkBBqXz58tq2bZu++OILbd26VTk5OXrwwQfVt2/fXIM7AAAAAEBRVOA7Sd3c3DRo0CANGjTInv0AAIqhGZ/9R85O/vLzTtfUUcvMbsdU9Xx36P2HpEi3HZJqm90OACAfBQpKn3zyyVXn9+/fv0DNAACKp5/XN1TiaT9VCDxX4oNSuPsR1WglZeQcMbsVAMBVFCgojRgxItf7rKwspaWlydnZWe7u7gQlAAAAAEVagUa9S05OzvVKTU3V3r171axZMwZzAAAAAFDkFSgo5SUiIkKvvPLKFWebAAAAAKCosVtQkqTSpUvr2LFj9iwJAAAAALdcge5RWrhwYa73hmEoISFBM2fOVNOmTe3SGAAAAACYpUBBqXv37rneWywWBQQEqHXr1po2bZo9+gIAAAAA0xQoKOXk5Ni7DwAAAAAoNOx6jxIAAAAAFAcFOqM0atSo6152+vTpBdkEAKAYaVZvh0qVKqey3ulmt2K6o+nB2rb7kO6oG6wws5sBAOSrQEHpjz/+0NatW3Xx4kVFRUVJkv766y+VLl1aMTEx1uUsFot9ugQAFGnPPjxXMdEhZrdRKPyeHKN+bxzSkjdjCEoAUIgVKCh17dpVXl5e+vjjj+Xn5yfp0kNoBw4cqDvvvFOjR4+2a5MAAAAAcCsV6B6ladOmafLkydaQJEl+fn566aWXGPUOAAAAQJFXoKB07tw5nThx4orpiYmJSklJuemmAAAAAMBMBbr0rkePHho4cKCmTZumO+64Q5K0ceNGjRkzRj179rRrgwCAoq/fM8/oXKqfgv1TtfmL98xux1Qdg1bqyJuSl/dKSbXNbgcAkI8CBaV33nlHTz75pPr166esrKxLhcqU0YMPPqipU6fatUEAQNF36oy3Ek97m91GoeBWOkPlykoZORlmtwIAuIoCBSV3d3e9/fbbmjp1qg4cOCDDMFStWjV5eHjYuz8AAAAAuOVu6oGzCQkJSkhIUGRkpDw8PGQYhr36AgAAAADTFCgonTp1Sm3atFFkZKQ6d+6shIQESdJDDz3E0OAAAAAAirwCBaUnnnhCTk5Oio+Pl7u7u3V6r169tGTJErs1BwAAAABmKNA9SkuXLtXPP/+s0NDQXNMjIiJ0+PBhuzQGAAAAAGYp0Bml8+fP5zqTdFlSUpJcXFxuuikAAAAAMFOBglLz5s31ySefWN9bLBbl5ORo6tSpatWqld2aAwAAAAAzFOjSu6lTp6ply5bavHmzLly4oLFjx2rXrl06ffq01q1bZ+8eAQAAAOCWKlBQqlGjhrZv365Zs2apdOnSOn/+vHr27KmhQ4cqJCTE3j0CAIq44X2+VVC5QLm7Zpndiun+OFNbs+dv0vD/1lYds5sBAOTrhoNSVlaW2rdvr3fffVcvvPCCI3oCABQznZptUkw0/5AmSYfSwvTh6k269z9hBCUAKMRu+B4lJycn7dy5UxaLxRH9AAAAAIDpCjSYQ//+/fXhhx/auxcAAAAAKBQKdI/ShQsX9MEHH2jZsmVq0KCBPDw8cs2fPn26XZoDABQPh44FycUpQGXK5Cgq/JTZ7ZjKq0yKalSQ3EulmN0KAOAqbigo/f333woPD9fOnTsVExMjSfrrr79yLcMleQAAW4+99IQST/upQuA5/bOsZP9jWtvAtbp7ipSRs1ZSE7PbAQDk44aCUkREhBISErRq1SpJUq9evfTGG28oKCjIIc0BAAAAgBlu6B4lwzByvV+8eLHOnz9v14YAAAAAwGwFGszhMtvgBAAAAADFwQ0FJYvFcsU9SNyTBAAAAKC4uaF7lAzD0IABA+Ti4iJJysjI0ODBg68Y9e67776zX4cAAAAAcIvdUFB64IEHcr3v16+fXZsBAAAAgMLghoLS7NmzHdUHAKAYy8nJ+f8/s5WQkGCXmsnJyXapAwBAXgr0wFkAAK7H5VCUnp4uSUpNPa9333vPLrV3xF/6mZbO6KsAAPszNShNnjxZ3333nfbs2SM3Nzc1adJEr776qqKioqzLGIahF154Qe+9956Sk5PVqFEjvfXWW6pZs6aJnQOAOeLj45WUlOSQ2v7+/qpYsaJda545c0aSVLr0pXtbnZ09VL/+I3apferiVkmbdSHzgl3qAQDwb6YGpTVr1mjo0KFq2LChLl68qGeffVbt27fX7t27rQNETJkyRdOnT9ecOXMUGRmpl156Se3atdPevXvl5eVlZvsAcEvFx8erevVopaenOaS+m5u79uyJs3tYkqRR9z6tGlVqqlSpHHl5hdilprNL0fx/wOLjrTTi1cX6dGIrtTK7GQBAvkwNSkuWLMn1fvbs2QoMDNSWLVvUvHlzGYahGTNm6Nlnn1XPnj0lSR9//LGCgoL0+eef69FHHzWjbQAwRVJSktLT09Sjx2cKCIi2a+2TJ+M0f34/JSUlOSQo+XicUYDfObvXLYoyctx09LR0wXAzuxUAwFUUqnuUzp49K0kqW7asJOngwYM6fvy42rdvb13GxcVFLVq00Pr16/MMSpmZmcrMzLS+P3eO/zEDKF4CAqIVEhJjdhsAABRrN/TAWUcyDEOjRo1Ss2bNVKtWLUnS8ePHJUlBQUG5lg0KCrLOszV58mT5+PhYX2FhYY5tHAAAAECxU2iC0uOPP67t27friy++uGKexWLJ9d4wjCumXTZu3DidPXvW+jpy5IhD+gUAXL/1u9rpq5WN9cO6+ma3YrpqHgf1RCcp1OWg2a0AAK6iUFx6N2zYMC1cuFC//PKLQkNDrdODg4MlXTqzFBLyfzf/JiYmXnGW6TIXFxe5uLg4tmEAwA1ZuvlenT1fTv4+59S16Raz2zFVbZ84NeonZeTESbrb7HYAAPkw9YySYRh6/PHH9d1332nlypWqXLlyrvmVK1dWcHCwli1bZp124cIFrVmzRk2aNLnV7QIAAAAoIUw9ozR06FB9/vnn+v777+Xl5WW978jHx0dubm6yWCwaOXKkJk2apIiICEVERGjSpElyd3dXnz59zGwdAAAAQDFmalCaNWuWJKlly5a5ps+ePVsDBgyQJI0dO1bp6ekaMmSI9YGzS5cu5RlKAAAAABzG1KBkGMY1l7FYLIqNjVVsbKzjGwIAAAAAFaJR7wAAAACgsCAoAQAAAICNQjE8OAAUJ/Hx8UpKSrJ73bi4OLvXdPQ2Dh7kWUEAgKKJoAQAdhQfH6/q1aOVnp7msG2kpqY4oGaCJIv69etn99qXXPueVAAAChOCEgDYUVJSktLT09Sjx2cKCIi2a+19+37SqlXjlZGRYde6kpSRcUaSoVatZioiorHd6m7d+oE2b56lct5H5euVo7Je5+1Wu6g6d9FLB49nKDjIS65mNwMAyBdBCQAcICAgWiEhMXatmZTk+Evv/Pyq2bVvL6/ykqSHO43THbfdbre6RdmKxDvVb9x3WvLmnepgdjMAgHwxmAMAAAAA2CAoAQAAAIANghIAAAAA2OAeJQCAw81bPUbfrg+Tj0eannvgO7PbMVWTcpu05CmptscmSbXNbgcAkA+CEgDA4Q4er61zaf7y9zlndiumC3I5qcp1pIyck2a3AgC4Ci69AwAAAAAbBCUAAAAAsEFQAgAAAAAbBCUAAAAAsMFgDgCAIi0lNUUJCQl2rZmcnGzXegCAooegBAAoki5mZ0iSNv2+Wf/s22zX2jviL/1MSz9v17oAgKKDoASgxIqPj1dSUpJda8bFxdm1HvKXczFLkhQQUFf16zaya+1TF7dK2qwLmRfsWhcAUHQQlACUSPHx8apePVrp6WkOqZ+amuKQuriSk7OHvLxC7FrT2cXLrvUAAEUPQQlAiZSUlKT09DT16PGZAgKi7VZ3376ftGrVeGVkZNitZnHQMGqJ/HyqytON/bI/tbKWr4/T3S0rq5rZzQAA8kVQAlCiBQREKyQkxm71kpK49C4vbet9rjtuu93sNgqFHeeiNXpunGreEU1QAoBCjOHBAQAAAMAGZ5QAAMgHQ48DQMlFUAIAwAZDjwMACEoAAIebPO8TnfvIX/4+5/T1i9PNbueaHDn0+IhOs/XtE1lKu7ha0h12rQ0AsB+CEgAA+XDE0OMWC7cHA0BRwN/WAAAAAGCDoAQAAAAANghKAAAAAGCDoAQAAAAANghKAAAAAGCDoAQAAAAANhgeHMBNi4+PV1JSkkNq+/v7q2LFig6pDQAAkB+CEoCbEh8fr+rVo5WenuaQ+m5u7tqzJ46wBAAAbimCEoCbkpSUpPT0NPXo8ZkCAqLtWvvkyTjNn99PSUlJBKUirleLqaoWXkvOZbLNbsV0n2+J1ve/btNT/aPV2uxmAAD5IigBsIuAgGiFhMSY3QYKqSohO3R7tJvZbRQKB5L8tHSHNCjDz+xWAABXwWAOAAAAAGCDoAQAAAAANrj0DgDgcH8n1FYpl6pyLpOt2yIOmd2Oqar6J6t9bSnANdnsVgAAV0FQAgA43JdrxujcYn/5+5zT1y9ON7sdU/WpH6chzaS0i3GSupndDgAgHwQlAIVeXFxckagJoOAc9TvJs9gAFBRBCUChlZqaIMmifv36OXAbKQ6rDeDaEpJSZZEc9nvu7u6muLg9hCUAN4ygBKDQysg4I8lQq1YzFRHR2K619+37SatWjVdGRoZd6wK4MWdSMmRImjm2lRrHRNi1dtzfJ9Xvmfk8iw1AgRCUABR6fn7V7P6MpqQkLr0DCpNqYX6KiQ4xuw0AsGJ4cAAAAACwQVACAAAAABsEJQAAAACwQVACAAAAABsEJQAAAACwwah3AACHG9e7v+647Xaz2ygUXvy5iaZ/vUrzXmyiXreZ3Q0AID+cUQIAAAAAGwQlAAAAALBBUAIAAAAAG9yjBABwuOV/9NGWv6vK0y1DD3RaY3Y7pmofdVAV+ko1fQ9Kam52OwCAfBCUAAAOt2lvR51L85e/z7kSH5QahSeoQ7SUdjHB7FYAAFdBUAIAWCUn71NCQoDd6qWkHLNbLQAAbiWCEgBAaWlJskhatWqYVq1yxBYMRxQFAMBhCEoAAGVmpsiQ9Hi3JqodWctudVf8vlZfrYmTQVACABQxBCUAgFWFcj6KDAuxW70dez3sVqu4MXJylJBgv/uUkpOT7VYLAEBQAgDg1jIunV27kJWld997z25ld8Rf+pmWft5uNQGgJCMoASVEfHy8kpKS7F43Li7O7jWB4uzyZYilSjmpfv2Bdqt76uJWSZt1IfOC3WoCQElGUAJKgPj4eFWvHq309DSHbSM1NcVhtYHiyGKxyMvLfpc5Ort42a0WAICgBJQISUlJSk9PU48enykgINqutfft+0mrVo1XRkaGXesCAACYiaAElCABAdEKCYmxa82kJC69w7WFB26Ti2sV+Xg47qxmUbH/pJdOnz2tEH9vycX+9VNSU+w6SITEQBEASiaCEgDA4f7T7AW1uL2d2W0UCp9srKKZC08rtl+kWtxuv7oXsy+d1d30+2b9s2+z/QqLgSIAlEwEJQAAioGci1mSpICAuqpft5FdazNQBICSiKAEAEAx4uTsYddBIqRbM1DEmTPJdr9k8ORJ+4/0CaDkICgBAADTpJ2/dN/aypWrFPfHKrvWTvj/t1bZO4ABKBkISgAAh5uz7A3NXh6msl7nNX3Yx2a3Y6phLfdqcFPJ3W2nDqto3bfliIEiTp26dNanfIUGiom272Azuw7sk1at0pkzZ+xaF0DJQFACADjcqZQwpaQF6nz6ObNbMV2AV4b83KVzGRk6XETGRrgVA0WULuVm90sGXV259A5AwRGUAADAVTlyoIiE9LWS4pSdk23XugBwswhKAADgujhkoAgnD7vWAwB7KWV2AwAAAABQ2BCUAAAAAMAGQQkAAAAAbBCUAAAAAMAGQQkAAAAAbBCUAAAAAMAGw4MDAByuRe3ZCgu5TW4uF8xuxXRLdpXXpr2H1eq2UNl5pG0AgB0RlAAADtcgYqFa3J5udhuFwvq/AzRz8WH5lAtWC4ISABRaXHoHAAAAADYISgAAAABgg6AEAHC4lLRyOpnsrVNnPc1uxXTerhdUoazkyf1aAFComRqUfvnlF3Xt2lXly5eXxWLRggULcs03DEOxsbEqX7683Nzc1LJlS+3atcucZgEABfbekg9034RRGvzaI2a3Yron28Xpnzelh5tsN7sVAMBVmBqUzp8/r7p162rmzJl5zp8yZYqmT5+umTNnatOmTQoODla7du2UkpJyizsFAAAAUJKYOupdp06d1KlTpzznGYahGTNm6Nlnn1XPnj0lSR9//LGCgoL0+eef69FHH72VrQIAAAAoQQrtPUoHDx7U8ePH1b59e+s0FxcXtWjRQuvXr893vczMTJ07dy7XCwAAAABuRKENSsePH5ckBQUF5ZoeFBRknZeXyZMny8fHx/oKCwtzaJ8AAAAAip9CG5Qus1gsud4bhnHFtH8bN26czp49a30dOXLE0S0CAAAAKGZMvUfpaoKDgyVdOrMUEvJ/jy5PTEy84izTv7m4uMjFxcXh/QEAAAAovgrtGaXKlSsrODhYy5Yts067cOGC1qxZoyZNmpjYGQAAAIDiztQzSqmpqdq/f7/1/cGDB7Vt2zaVLVtWFStW1MiRIzVp0iRFREQoIiJCkyZNkru7u/r06WNi1wAAAACKO1OD0ubNm9WqVSvr+1GjRkmSHnjgAc2ZM0djx45Venq6hgwZouTkZDVq1EhLly6Vl5eXWS0DAAAAKAFMDUotW7aUYRj5zrdYLIqNjVVsbOytawoAYHcPtBmhmFrNVLpUjtmtmG7m6ijNW71Tj3WpqVo1ze4GAJCfQjuYAwCg+PD3iVflkJNmt1EoJKa4avdR6dR5N7NbAQBcRaEdzAEAAAAAzMIZJaCQiY+PV1JSkl1rxsXF2bUeAABAcUdQAgqR+Ph4Va8erfT0NIfUT01NcUhd4Fq2H2yn1IsxcnHOUtsGO8xux1T1K57Sgy2lWlyKCACFGkEJKESSkpKUnp6mHj0+U0BAtN3q7tv3k1atGq+MjAy71QRuxLI/hihlXaD8fc6V+KDUre4/GtBYOpdxWFvPm90NACA/BCWgEAoIiFZISIzd6iUlcekdAADAjWAwBwAAAACwQVACAAAAABsEJQAAAACwQVACAAAAABsEJQAAAACwQVACAAAAABsEJQAAAACwwXOUAAAO5+l6Wi5OrirrnWp2K6ZLyXDS+YwslSrtZHYrAICrICgBABzu0c4PqsXt7cxuo1CYuqyGZi7crNh+ddXidrO7AQDkh0vvAAAAAMAGQQkAAAAAbHDpHQAAKNYOHjyorVu32rWmv7+/KlasaNeaAAoXghIAwOF+2DhGq3ZGyMs9XaN7/2h2O6bq1eCQmodLVUIOKEXct+VIWRfPS5LGjx+v8ePH27W2m5u79uyJIywBxRhBCQDgcH8da6KUtED5+5wzuxXT1Qw5q2ZVpXMZydp63uxuirfsixckSQ0aPKWYmPvsVvfkyTjNn99PSUlJBCWgGCMoAQCAYs3LK1QhITFmtwGgiGEwBwAAAACwwRklAACAAoiLi3NIXQaKAAoHghIAAMANSE1NkGRRv379HFKfgSKAwoGgBAAAcAMyMs5IMtSq1UxFRDS2a20GigAKD4ISAAAo1lJSjighwX7PUUpOPihJ8vOrxiARQDFGUAIAAMXSmdQMWSRt3jxFmzdPsWtti6S0tES71gRQuBCUAABAsXQ+I0uGpEc611H9mnfYre6Ov7Zq5veblZl51m41ARQ+BCUAgMPVrrRcvj7V5eWebnYrptsSX1aHEk6oRiV/ycfsbkqG8mU9FRkWYrd6p0952a0WgMKLoAQAcLj29d9Si9vbmd1GofD9n2GaufCEYvuFq8XtZncDAMgPD5wFAAAAABsEJQAAAACwQVACAAAAABvcowQAcLg3F36uKV8Hq5xPij55bqbZ7ZjquU479fLd0oXsP7Qzk/u2kLe4uDiH1PX39+dBtsB1IigBABzuwkU3pWW6yD0z0+xWTOdcJlvebtK5jGyJ3QEbqakJkizq16+fQ+q7ublrz544whJwHQhKAAAAhURGxhlJhlq1mqmIiMZ2rX3yZJzmz++npKQkghJwHQhKAAAAhYyfXzWFhMSY3QZQojGYAwAAAADY4IwSABQxycn7lJAQYNeaKSnH7FoPAICijqAEAEVEWlqSLJJWrRqmVascs42LF9MdUxgAgCKGoAQARURmZooMSY93a6LakbXsWnvF72v11Zo4ZedcsGtdAACKKoISABQxFcr5KDIsxK41d+z1sGs9AACKOgZzAAAAAAAbnFECADhcl9unKrJyQ7k4XTS7FdN9ubmSVvyxXz2aVlH5ymZ3AwDID0EJAOBwUaHr1aIel/dJ0q4EX33zu1QrsixBCQAKMS69AwAAAAAbBCUAAAAAsEFQAgA43LFTUdp1MFR74+07Wl9RFOZ3XndUk0K8U81uBQBwFdyjBABwuC/WvKL3FgfK3+ecvn5xutntmOrhZvs1tr10LmOPtp43uxsAQH44owQAAAAANghKAAAAAGCDoAQAAAAANghKAAAAAGCDoAQAAAAANghKAAAAAGCD4cGBAoiPj1dSUpLd68bFxdm9JgDAMVJSjighYatdayYnH7RrPQAFR1ACblB8fLyqV49Wenqaw7aRmprisNoAgJtzJjVDFkmbN0/R5s1T7F7fIiktLdHudQHcGIIScIOSkpKUnp6mHj0+U0BAtF1r79v3k1atGq+MjAy71gUA2M/5jCwZkh7pXEf1a95h19o7/tqqmd9vVmbmWbvWBXDjCEpAAQUERCskJMauNZOSuPQOxdPjXfuqaUxrWSyG2a2Y7uXFtfTe4j807r56alzf7G5wM8qX9VRkWIhda54+5WXXegAKjqAEAHA4F6c0ebhlmt1GoZB5sbRS0qUL2aXNbgUAcBWMegcAAAAANjijBFM5avQ4ScrMzJSLi4vd6zIyHQAAQPFHUIJpHD96nEWS4+6HYGQ64Pqt391LB5PqyN01U/e13mB2O6ZqFXlc5XpKd4QfM7sVAMBVEJRgmlsxelyrVjMVEdHYIbUZmQ64fhv29FbK1kD5+5wjKEWdUM960rmMY9p63uxuAAD5ISjBdI4cPc7Prxoj0wEAAOCGMZgDAAAAANggKAEAAACADYISAAAAANjgHiUAcIDk5H1KSAiwa82UFEZJA3DzHPGYC0c9kkOS/P39VbFiRYfUBq6GoAQAdpSWliSLpFWrhmnVKsds4+LFdMcUBlCspaYmSLKoX79+DqjuuEdyuLm5a8+eOMISbjmCEgDYUWZmigxJj3drotqRtexae8Xva/XVmjhl51ywa10AJUNGxhlJht0fneHIR3KcPBmn+fP7KSkpiaCEW46gBAAOUKGcjyLDQuxac8deD7vWA1Ay2fvRGY58JMdljrhcUOKyPlwdQQkA4HAhfn8pLDBTvp5pZrdiun+S3bXn6FmV9fbg/8LANTj2ckEu68PV8Vc0AMDh+rR6Si1ub2d2G4XCe79GaObCzYrtF60Wt5vdDVC4OepyQYnL+nBtBCUAAAAUao68rA/ID89RAgAAAAAbnFEqRuLj45WUlOSQ2tzsCAAAgJKEoFRMxMfHq3r1aKWnO+ZGaW52BHAzPl/1qr76tZJ8PdP08iNfmN2OqR5ptk99YqSy3nFKEPdtAUBhRVC6xRx11icuLk7p6Wnq0eMzBQRE27X25Zsd165dq+ho+9V21FCfAAqfhORIpaQFyt/nnNmtmC7UL01+FaRzGeeVcN7sbgA4iqO+83GVz61DULqFHH3WR5Lc3ELtfrOjo4fmTE1NcUhdAAAAMzjyOx9X+dw6BKVbKCkpyWFnfS4/FTsjI8OudSXHP8nbET0DAACYxVHf+RjS/NYiKJkgICDa7md9Lj8V25Ec9SRvAACQW0rKESUkbLVrzeTkg3atV1w44laAyzUd8Z3v3/Xtjcv6ciMoAQAAFBJnUjNkkbR58xRt3jzF7vUtktLSEu1etyhy9K0Fl7Zh39sLHN0zl/XlViSC0ttvv62pU6cqISFBNWvW1IwZM3TnnXea3RYAAIBdnc/IkiHpkc51VL/mHXatveOvrZr5/WZlZp61a92iylG3FkiOu73AkT1zWd+VCn1Q+vLLLzVy5Ei9/fbbatq0qd5991116tRJu3fv5iACAIBiqXxZT0WGhdi15ulTXnatV1zY+9YCyfG3FziiZ1yplNkNXMv06dP14IMP6qGHHlJ0dLRmzJihsLAwzZo1y+zWAAAAABRThfqM0oULF7RlyxY9/fTTuaa3b99e69evz3OdzMxMZWZmWt+fPXvp9PK5c+Y/uyM1NVWSFBf3nRIStti19j///CZJOnhwibKy4otE7aLYc1GtXRR7Lqq1jx7dLEnaeyReKmXfv2IPn7h0X8H+o8flst2+N3k7qvbluhcvnpXkqgtZp7XRTvWL4v6QpPBa2Sot6XRatjbusP++Lmr7oyjWLoo9S9LeI0clSUePrtfWrS52q1sU/64uqrUd2fOZM4clSd999522bLHv91RJio6OVp06dexe90ZdzgSGYVx7YaMQO3r0qCHJWLduXa7pL7/8shEZGZnnOhMmTDAk8eLFixcvXrx48eLFi1eeryNHjlwzixTqM0qXWSyWXO8Nw7hi2mXjxo3TqFGjrO9zcnJ0+vRplStXLt914Fjnzp1TWFiYjhw5Im9vb7PbgYNxvEsOjnXJwbEuOTjWJUdJPdaGYSglJUXly5e/5rKFOij5+/urdOnSOn78eK7piYmJCgoKynMdFxcXubjkPp3s6+vrqBZxA7y9vUvUL2JJx/EuOTjWJQfHuuTgWJccJfFY+/j4XNdyhXowB2dnZ9WvX1/Lli3LNX3ZsmVq0qSJSV0BAAAAKO4K9RklSRo1apTuv/9+NWjQQI0bN9Z7772n+Ph4DR482OzWAAAAABRThT4o9erVS6dOndLEiROVkJCgWrVq6aefflKlSpXMbg3XycXFRRMmTLjikkgUTxzvkoNjXXJwrEsOjnXJwbG+NothXM/YeAAAAABQchTqe5QAAAAAwAwEJQAAAACwQVACAAAAABsEJQAAAACwQVDCdfnll1/UtWtXlS9fXhaLRQsWLMg132Kx5PmaOnWqdZnMzEwNGzZM/v7+8vDw0N13361//vknV53k5GTdf//98vHxkY+Pj+6//36dOXPmFnxCXHazx/r06dMaNmyYoqKi5O7urooVK2r48OE6e/Zsrjoc68LBHr/blxmGoU6dOuVZh+NtPnsd6w0bNqh169by8PCQr6+vWrZsqfT0dOt8jrX57HGsjx8/rvvvv1/BwcHy8PBQTEyMvvnmm1x1ONbmu9axTk1N1eOPP67Q0FC5ubkpOjpas2bNyrUM38/yR1DCdTl//rzq1q2rmTNn5jk/ISEh1+ujjz6SxWLRf/7zH+syI0eO1Pz58zVv3jz9+uuvSk1NVZcuXZSdnW1dpk+fPtq2bZuWLFmiJUuWaNu2bbr//vsd/vnwf272WB87dkzHjh3Ta6+9ph07dmjOnDlasmSJHnzwwVx1ONaFgz1+ty+bMWOGLBZLnnU43uazx7HesGGDOnbsqPbt2+v333/Xpk2b9Pjjj6tUqf/7OsGxNp89jvX999+vvXv3auHChdqxY4d69uypXr166Y8//rAuw7E237WO9RNPPKElS5bos88+U1xcnJ544gkNGzZM33//vXUZvp9dhQHcIEnG/Pnzr7pMt27djNatW1vfnzlzxnBycjLmzZtnnXb06FGjVKlSxpIlSwzDMIzdu3cbkoyNGzdal9mwYYMhydizZ499PwSuS0GOdV6++uorw9nZ2cjKyjIMg2NdWN3M8d62bZsRGhpqJCQkXFGH4134FPRYN2rUyHjuuefyXYdjXfgU9Fh7eHgYn3zySa5pZcuWNT744APDMDjWhVFex7pmzZrGxIkTc02LiYmx/h7z/ezqOKMEuztx4oQWLVqU6wzCli1blJWVpfbt21unlS9fXrVq1dL69eslXfqXSh8fHzVq1Mi6zB133CEfHx/rMihc8jrWeTl79qy8vb1VpsylZ1xzrIum/I53Wlqa/vvf/2rmzJkKDg6+Yj2Od9GT17FOTEzUb7/9psDAQDVp0kRBQUFq0aKFfv31V+syHOuiJ7/f62bNmunLL7/U6dOnlZOTo3nz5ikzM1MtW7aUxLEuKpo1a6aFCxfq6NGjMgxDq1at0l9//aUOHTpI4vvZtRCUYHcff/yxvLy81LNnT+u048ePy9nZWX5+frmWDQoK0vHjx63LBAYGXlEvMDDQugwKl7yOta1Tp07pxRdf1KOPPmqdxrEumvI73k888YSaNGmibt265bkex7voyetY//3335Kk2NhYPfzww1qyZIliYmLUpk0b7du3TxLHuijK7/f6yy+/1MWLF1WuXDm5uLjo0Ucf1fz581W1alVJHOui4o033lCNGjUUGhoqZ2dndezYUW+//baaNWsmie9n11LG7AZQ/Hz00Ufq27evXF1dr7msYRi57mnI6/4G22VQeFzrWJ87d0533XWXatSooQkTJuSax7EuevI63gsXLtTKlStz3beQF4530ZLXsc7JyZEkPfrooxo4cKAkqV69elqxYoU++ugjTZ48WRLHuqjJ7+/x5557TsnJyVq+fLn8/f21YMEC3XvvvVq7dq1q164tiWNdFLzxxhvauHGjFi5cqEqVKumXX37RkCFDFBISorZt2+a7Ht/PLiEowa7Wrl2rvXv36ssvv8w1PTg4WBcuXFBycnKuf7VITExUkyZNrMucOHHiiponT55UUFCQYxvHDcvvWF+WkpKijh07ytPTU/Pnz5eTk5N1Hse66MnveK9cuVIHDhyQr69vrun/+c9/dOedd2r16tUc7yImv2MdEhIiSapRo0au6dHR0YqPj5fE73ZRk9+xPnDggGbOnKmdO3eqZs2akqS6detq7dq1euutt/TOO+9wrIuA9PR0PfPMM5o/f77uuusuSVKdOnW0bds2vfbaa2rbti3fz66BS+9gVx9++KHq16+vunXr5ppev359OTk5admyZdZpCQkJ2rlzp/UXsXHjxjp79qx+//136zK//fabzp49a10GhUd+x1q6dCapffv2cnZ21sKFC6/4l0qOddGT3/F++umntX37dm3bts36kqTXX39ds2fPlsTxLmryO9bh4eEqX7689u7dm2v6X3/9pUqVKkniWBc1+R3rtLQ0Sco1mqEklS5d2npmkWNd+GVlZSkrK+uqx5HvZ9dg2jASKFJSUlKMP/74w/jjjz8MScb06dONP/74wzh8+LB1mbNnzxru7u7GrFmz8qwxePBgIzQ01Fi+fLmxdetWo3Xr1kbdunWNixcvWpfp2LGjUadOHWPDhg3Ghg0bjNq1axtdunRx+OfD/7nZY33u3DmjUaNGRu3atY39+/cbCQkJ1hfHuvCxx++2LeUx8hLH23z2ONavv/664e3tbXz99dfGvn37jOeee85wdXU19u/fb12GY22+mz3WFy5cMKpVq2bceeedxm+//Wbs37/feO211wyLxWIsWrTIuhzH2nzXOtYtWrQwatasaaxatcr4+++/jdmzZxuurq7G22+/ba3B97P8EZRwXVatWmVIuuL1wAMPWJd59913DTc3N+PMmTN51khPTzcef/xxo2zZsoabm5vRpUsXIz4+Ptcyp06dMvr27Wt4eXkZXl5eRt++fY3k5GQHfjLYutljnd/6koyDBw9al+NYFw72+N22lVdQ4nibz17HevLkyUZoaKjh7u5uNG7c2Fi7dm2u+Rxr89njWP/1119Gz549jcDAQMPd3d2oU6fOFcOFc6zNd61jnZCQYAwYMMAoX7684erqakRFRRnTpk0zcnJyrDX4fpY/i2EYhmPOVQEAAABA0cQ9SgAAAABgg6AEAAAAADYISgAAAABgg6AEAAAAADYISgAAAABgg6AEAAAAADYISgAAAABgg6AEAAAAADYISgCAW8pisWjBggVmtyGpcPUCAChcCEoAALs5fvy4hg0bpipVqsjFxUVhYWHq2rWrVqxYYdftDBgwQN27d7drzfw4+jPNmTNHvr6+dqkFALCfMmY3AAAoHg4dOqSmTZvK19dXU6ZMUZ06dZSVlaWff/5ZQ4cO1Z49e8xu8YYVx88EALg+nFECANjFkCFDZLFY9Pvvv+uee+5RZGSkatasqVGjRmnjxo35rnf06FH16tVLfn5+KleunLp166ZDhw7d0LZbtmyp4cOHa+zYsSpbtqyCg4MVGxuba5l9+/apefPmcnV1VY0aNbRs2TK7fab4+Hh169ZNnp6e8vb21n333acTJ05Y5//5559q1aqVvLy85O3trfr162vz5s1avXq1Bg4cqLNnz8pischisVzRNwDAHAQlAMBNO336tJYsWaKhQ4fKw8Pjivn5XVqWlpamVq1aydPTU7/88ot+/fVXeXp6qmPHjrpw4cIN9fDxxx/Lw8NDv/32m6ZMmaKJEydaw1BOTo569uyp0qVLa+PGjXrnnXf01FNP2eUzGYah7t276/Tp01qzZo2WLVumAwcOqFevXtZl+/btq9DQUG3atElbtmzR008/LScnJzVp0kQzZsyQt7e3EhISlJCQoCeffPKGPjcAwDG49A4AcNP2798vwzBUvXr1G1pv3rx5KlWqlD744ANZLBZJ0uzZs+Xr66vVq1erffv2112rTp06mjBhgiQpIiJCM2fO1IoVK9SuXTstX75ccXFxOnTokEJDQyVJkyZNUqdOnW76My1fvlzbt2/XwYMHFRYWJkn69NNPVbNmTW3atEkNGzZUfHy8xowZY60VERFhXd/Hx0cWi0XBwcHX/VkBAI7HGSUAwE0zDEOSrGHnem3ZskX79++Xl5eXPD095enpqbJlyyojI0MHDhzQ2rVrrdM9PT01d+7cfGvVqVMn1/uQkBAlJiZKkuLi4lSxYkVrSJKkxo0b2+UzxcXFKSwszBqSJKlGjRry9fVVXFycJGnUqFF66KGH1LZtW73yyis6cODAVWsCAMzHGSUAwE2LiIiQxWJRXFzcDY1Gl5OTo/r16+cZgAICAuTs7Kxt27ZZpwUFBeVby8nJKdd7i8WinJwcSf8XemznX831fibDMPKs9e/psbGx6tOnjxYtWqTFixdrwoQJmjdvnnr06HHVHgAA5uGMEgDgppUtW1YdOnTQW2+9pfPnz18x/8yZM3muFxMTo3379ikwMFDVqlXL9fLx8ZGbm1uuaV5eXgXqr0aNGoqPj9exY8es0zZs2GCXz3S59pEjR6zzdu/erbNnzyo6Oto6LTIyUk888YSWLl2qnj17avbs2ZIkZ2dnZWdnF+hzAQAch6AEALCLt99+W9nZ2br99tv17bffat++fYqLi9Mbb7yR72Vuffv2lb+/v7p166a1a9fq4MGDWrNmjUaMGKF//vnHbr21bdtWUVFR6t+/v/7880+tXbtWzz77rF0+U9u2bVWnTh317dtXW7du1e+//67+/furRYsWatCggdLT0/X4449r9erVOnz4sNatW6dNmzZZQ1R4eLhSU1O1YsUKJSUlKS0tzW6fGwBQcAQlAIBdVK5cWVu3blWrVq00evRo1apVS+3atdOKFSs0a9asPNdxd3fXL7/8oooVK6pnz56Kjo7WoEGDlJ6eLm9vb7v1VqpUKc2fP1+ZmZm6/fbb9dBDD+nll1+2y2eyWCxasGCB/Pz81Lx5c7Vt21ZVqlTRl19+KUkqXbq0Tp06pf79+ysyMlL33XefOnXqpBdeeEGS1KRJEw0ePFi9evVSQECApkyZYrfPDQAoOIuR14XbAAAAAFCCcUYJAAAAAGwQlAAAAADABkEJAAAAAGwQlAAAAADABkEJAAAAAGwQlAAAAADABkEJAAAAAGwQlADgOqxevVoWi0Vnzpyxe+3LDyzNT8uWLTVy5Ejr+/DwcM2YMcPufdjDoUOHZLFYtG3bNrNbKTZiY2N12223md0GAJQ4BCUARdqAAQNksVhksVjk5OSkoKAgtWvXTh999JFycnLMbs8hNm3apEceecTsNuAAeYXmJ598UitWrDCnIQAowQhKAIq8jh07KiEhQYcOHdLixYvVqlUrjRgxQl26dNHFixfNbs/uAgIC5O7ubnYbRUpWVpbZLRSYp6enypUrZ3YbAFDiEJQAFHkuLi4KDg5WhQoVFBMTo2eeeUbff/+9Fi9erDlz5liXmz59umrXri0PDw+FhYVpyJAhSk1Ntc4/fPiwunbtKj8/P3l4eKhmzZr66aefcm1ry5YtatCggdzd3dWkSRPt3bs31/wffvhB9evXl6urq6pUqaIXXnghV1jbt2+fmjdvLldXV9WoUUPLli274c9re+mdxWLRBx98oB49esjd3V0RERFauHBhrnV2796tzp07y9PTU0FBQbr//vuVlJRknf/NN9+odu3acnNzU7ly5dS2bVudP38+z+0nJyerb9++CggIkJubmyIiIjR79uxcy/z9999q1aqV3N3dVbduXW3YsCHX/G+//VY1a9aUi4uLwsPDNW3aNOu8N998U7Vr17a+X7BggSwWi9566y3rtA4dOmjcuHF59nf58r+vvvpKLVu2lKurqz777DNJ0uzZsxUdHS1XV1dVr15db7/9tnW9Cxcu6PHHH1dISIhcXV0VHh6uyZMn59rPs2bNUqdOneTm5qbKlSvr66+/zrXto0ePqlevXvLz81O5cuXUrVs3HTp0KNcyH330kfWzh4SE6PHHH5d06bhKUo8ePWSxWKzv/33p3c8//yxXV9crLgEdPny4WrRoYX2/fv16NW/eXG5ubgoLC9Pw4cPzPZ7/3sZHH32kihUrytPTU4899piys7M1ZcoUBQcHKzAwUC+//HKu9c6ePatHHnlEgYGB8vb2VuvWrfXnn39a5x84cEDdunVTUFCQPD091bBhQy1fvjxXjfDwcE2aNEmDBg2Sl5eXKlasqPfee++6jwsAOIwBAEXYAw88YHTr1i3PeXXr1jU6depkff/6668bK1euNP7++29jxYoVRlRUlPHYY49Z5991111Gu3btjO3btxsHDhwwfvjhB2PNmjWGYRjGqlWrDElGo0aNjNWrVxu7du0y7rzzTqNJkybW9ZcsWWJ4e3sbc+bMMQ4cOGAsXbrUCA8PN2JjYw3DMIzs7GyjVq1aRsuWLY0//vjDWLNmjVGvXj1DkjF//vx8P2OLFi2MESNGWN9XqlTJeP31163vJRmhoaHG559/buzbt88YPny44enpaZw6dcowDMM4duyY4e/vb4wbN86Ii4sztm7darRr185o1aqVdX6ZMmWM6dOnGwcPHjS2b99uvPXWW0ZKSkqe/QwdOtS47bbbjE2bNhkHDx40li1bZixcuNAwDMM4ePCgIcmoXr268eOPPxp79+417rnnHqNSpUpGVlaWYRiGsXnzZqNUqVLGxIkTjb179xqzZ8823NzcjNmzZxuGYRjbt283LBaLcfLkScMwDGPkyJGGv7+/ce+99xqGYRhZWVmGp6ensXjx4jz7u9xDeHi48e233xp///23cfToUeO9994zQkJCrNO+/fZbo2zZssacOXMMwzCMqVOnGmFhYcYvv/xiHDp0yFi7dq3x+eef59rP5cqVM95//31j7969xnPPPWeULl3a2L17t2EYhnH+/HkjIiLCGDRokLF9+3Zj9+7dRp8+fYyoqCgjMzPTMAzDePvttw1XV1djxowZxt69e43ff//deiwTExMNScbs2bONhIQEIzEx0TAMw5gwYYJRt25dwzAM4+LFi0ZQUJDxwQcfWPu6PO3dd9+17j9PT0/j9ddfN/766y9j3bp1Rr169YwBAwbkub8ub8PT09O45557jF27dhkLFy40nJ2djQ4dOhjDhg0z9uzZY3z00UeGJGPDhg2GYRhGTk6O0bRpU6Nr167Gpk2bjL/++ssYPXq0Ua5cOeufvW3bthnvvPOOsX37duOvv/4ynn32WcPV1dU4fPiwdduVKlUyypYta7z11lvGvn37jMmTJxulSpUy4uLiruu4AICjEJQAFGlXC0q9evUyoqOj8133q6++MsqVK2d9X7t2bWuosXU5KC1fvtw6bdGiRYYkIz093TAMw7jzzjuNSZMm5Vrv008/NUJCQgzDMIyff/7ZKF26tHHkyBHr/MWLF9slKD333HPW96mpqYbFYrEGifHjxxvt27fPVfPIkSOGJGPv3r3Gli1bDEnGoUOH8u3h37p27WoMHDgwz3mXQ8q/v8jv2rXLkGT94tunTx+jXbt2udYbM2aMUaNGDcMwLn0B9/f3N7755hvDMAzjtttuMyZPnmwEBgYahmEY69evN8qUKZNvkLvcw4wZM3JNDwsLu+IL9osvvmg0btzYMAzDGDZsmNG6dWsjJycnz7qSjMGDB+ea1qhRI2vY/vDDD42oqKhc62dmZhpubm7Gzz//bBiGYZQvX9549tln86x/eRu2fxb+HZQMwzCGDx9utG7d2vr+559/NpydnY3Tp08bhmEY999/v/HII4/kqrF27VqjVKlS1j+rtiZMmGC4u7sb586ds07r0KGDER4ebmRnZ1unRUVFGZMnTzYMwzBWrFhheHt7GxkZGblqVa1a1Rra8lKjRg3jzTfftL6vVKmS0a9fP+v7nJwcIzAw0Jg1a5ZhGNc+LgDgKFx6B6DYMgxDFovF+n7VqlVq166dKlSoIC8vL/Xv31+nTp2yXpI0fPhwvfTSS2ratKkmTJig7du3X1GzTp061v8OCQmRJCUmJkq6dFnexIkT5enpaX09/PDDSkhIUFpamuLi4lSxYkWFhoZaazRu3Ngun/XffXl4eMjLyytXX6tWrcrVV/Xq1SVdujSqbt26atOmjWrXrq17771X77//vpKTk/Pd1mOPPaZ58+bptttu09ixY7V+/fqr9mO7n+Li4tS0adNcyzdt2lT79u1Tdna2LBaLmjdvrtWrV+vMmTPatWuXBg8erOzsbMXFxWn16tWKiYmRp6fnVfdJgwYNrP998uRJHTlyRA8++GCu/fDSSy/pwIEDki4NDLJt2zZFRUVp+PDhWrp06RU1bY9X48aNFRcXJ+nSft6/f7+8vLys9cuWLauMjAwdOHBAiYmJOnbsmNq0aXPVvq+lb9++Wr16tY4dOyZJmjt3rjp37iw/Pz9rH3PmzMn1OTt06KCcnBwdPHgw37rh4eHy8vKyvg8KClKNGjVUqlSpXNP+/ecqNTVV5cqVy7WtgwcPWvfp+fPnNXbsWNWoUUO+vr7y9PTUnj17FB8fn2vb//7zYrFYFBwcbN3O9RwXAHCEMmY3AACOEhcXp8qVK0u6dP9R586dNXjwYL344osqW7asfv31Vz344IPWG/0feughdejQQYsWLdLSpUs1efJkTZs2TcOGDbPWdHJysv735RB2eXS9nJwcvfDCC+rZs+cVvbi6usowjCum/zvI3Yx/93W57r/76tq1q1599dUr1gsJCVHp0qW1bNkyrV+/XkuXLtWbb76pZ599Vr/99pt1//1bp06ddPjwYS1atEjLly9XmzZtNHToUL322mt59mO7n2wD7OVp/9ayZUu99957Wrt2rerWrStfX181b95ca9as0erVq9WyZctr7hMPDw/rf1/e9vvvv69GjRrlWq506dKSpJiYGB08eFCLFy/W8uXLdd9996lt27b65ptvrrqdf3+++vXra+7cuVcsExAQkCtw3Izbb79dVatW1bx58/TYY49p/vz5ue4Ry8nJ0aOPPqrhw4dfsW7FihXzrZvXn6Fr/bkKCQnR6tWrr6jl6+srSRozZox+/vlnvfbaa6pWrZrc3Nx0zz336MKFC9fc9uXtFPS4AMDNIigBKJZWrlypHTt26IknnpAkbd68WRcvXtS0adOsX1i/+uqrK9YLCwvT4MGDNXjwYI0bN07vv/9+rqB0NTExMdq7d6+qVauW5/waNWooPj5ex44dU/ny5SXpikEOHCEmJkbffvutwsPDVaZM3n/tWywWNW3aVE2bNtXzzz+vSpUqaf78+Ro1alSeywcEBGjAgAEaMGCA7rzzTo0ZMyZXULqaGjVq6Ndff801bf369YqMjLSGlpYtW2rEiBH65ptvrKGoRYsWWr58udavX68RI0Zc56e/JCgoSBUqVNDff/+tvn375ruct7e3evXqpV69eumee+5Rx44ddfr0aZUtW1aStHHjRvXv39+6/MaNG1WvXj1Jl/bzl19+aR3YIC/h4eFasWKFWrVqled8JycnZWdnX/Pz9OnTR3PnzlVoaKhKlSqlu+66yzovJiZGu3btyvfPob3ExMTo+PHjKlOmjHXgCVtr167VgAED1KNHD0lSamrqFYNbXI9rHRcAcASCEoAiLzMzU8ePH1d2drZOnDihJUuWaPLkyerSpYv1S23VqlV18eJFvfnmm+ratavWrVund955J1edkSNHqlOnToqMjFRycrJWrlyp6Ojo6+7j+eefV5cuXRQWFqZ7771XpUqV0vbt27Vjxw699NJLatu2raKiotS/f39NmzZN586d07PPPmvXfZGXoUOH6v3339d///tfjRkzRv7+/tq/f7/mzZun999/X5s3b9aKFSvUvn17BQYG6rffftPJkyfz/ezPP/+86tevr5o1ayozM1M//vjjDe2n0aNHq2HDhnrxxRfVq1cvbdiwQTNnzsw1Al2tWrVUrlw5zZ07V99//72kS+Fp9OjRkqRmzZrd8H6IjY3V8OHD5e3trU6dOikzM1ObN29WcnKyRo0apddff10hISG67bbbVKpUKX399dcKDg62nh2RpK+//loNGjRQs2bNNHfuXP3+++/68MMPJV26JG7q1Knq1q2bJk6cqNDQUMXHx+u7777TmDFjFBoaqtjYWA0ePFiBgYHq1KmTUlJStG7dOmsYvxykmjZtKhcXF+vldLb69u2rF154QS+//LLuueceubq6Wuc99dRTuuOOOzR06FA9/PDD8vDwUFxcnJYtW6Y333zzhvdbftq2bavGjRure/fuevXVVxUVFaVjx47pp59+Uvfu3dWgQQNVq1ZN3333nbp27SqLxaLx48ff8PPNrue4AIAjcI8SgCJvyZIlCgkJUXh4uDp27KhVq1bpjTfe0Pfff289Q3Hbbbdp+vTpevXVV1WrVi3NnTv3iiGGs7OzNXToUEVHR6tjx46KiorK9eX9Wjp06KAff/xRy5YtU8OGDXXHHXdo+vTpqlSpkiSpVKlSmj9/vjIzM3X77bfroYceumK4ZUcoX7681q1bp+zsbHXo0EG1atXSiBEj5OPjo1KlSsnb21u//PKLOnfurMjISD333HOaNm2aOnXqlGc9Z2dnjRs3TnXq1FHz5s1VunRpzZs377r7iYmJ0VdffaV58+apVq1aev755zVx4kQNGDDAuozFYrEOd33nnXdKunQfi4+Pj+rVq5fvGZureeihh/TBBx9ozpw5ql27tlq0aKE5c+ZYLy/09PTUq6++qgYNGqhhw4Y6dOiQfvrpp1yXzL3wwguaN2+e6tSpo48//lhz585VjRo1JEnu7u765ZdfVLFiRfXs2VPR0dEaNGiQ0tPTrf0+8MADmjFjht5++23VrFlTXbp00b59+6z1p02bpmXLliksLMx6piovERERatiwobZv337FGbI6depozZo12rdvn+68807Vq1dP48ePt94rZi8Wi0U//fSTmjdvrkGDBikyMlK9e/fWoUOHFBQUJOlSyPHz81OTJk3UtWtXdejQQTExMTe0nes5LgDgCBYjr4vmAQBALhaLRfPnz1f37t3NbgUAcAvwzzEAAAAAYIOgBAAAAAA2GMwBAIDrwJXqAFCycEYJAAAAAGwQlAAAAADABkEJAAAAAGwQlAAAAADABkEJAAAAAGwQlAAAAADABkEJAAAAAGwQlAAAAADABkEJAAAAAGz8P8dtf1+JZSk/AAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 1000x600 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "np.random.seed(2579)\n",
    "\n",
    "# Number of samples to generate\n",
    "n_samples = 1000\n",
    "\n",
    "# Current cost between Cleveland and Indianapolis\n",
    "mean_cost = 2.3\n",
    "std_sev = 0.2\n",
    "# Generate random costs using normal distribution\n",
    "normal_costs = np.random.normal(loc=mean_cost, scale=std_sev, size=n_samples)\n",
    "\n",
    "# Extract the shipment amount from Cleveland to Indianapolis for each solution\n",
    "shipment_amount = sol_sim_comp.loc[('Cleveland', 'Indianapolis'), 'shipment']\n",
    "shipment_amount2 = sol_sim_comp.loc[('Cleveland', 'Indianapolis'), 'shipment2']\n",
    "\n",
    "# Find the total costs for all other shipments for each of the two solution since we are assuming they are fixed\n",
    "obj1 = sum(sol_sim_comp.cost * sol_sim_comp.shipment) - shipment_amount*sol_sim_comp.loc[('Cleveland', 'Indianapolis'), 'cost']\n",
    "obj2 = sum(sol_sim_comp.cost * sol_sim_comp.shipment2) - shipment_amount2*sol_sim_comp.loc[('Cleveland', 'Indianapolis'), 'cost']\n",
    "\n",
    "# Create a DataFrame to hold total costs including the new simulated costs \n",
    "total_costs_df = pd.DataFrame({\n",
    "    'Total_Cost_Sol1': obj1 + normal_costs * shipment_amount,\n",
    "    'Total_Cost_Sol2': obj2 + normal_costs * shipment_amount2\n",
    "})\n",
    "\n",
    "# Determine the common range for both datasets\n",
    "min_cost = min(total_costs_df['Total_Cost_Sol1'].min(), total_costs_df['Total_Cost_Sol2'].min())\n",
    "max_cost = max(total_costs_df['Total_Cost_Sol1'].max(), total_costs_df['Total_Cost_Sol2'].max())\n",
    "\n",
    "### This is code to create a visualization of each solution's variability\n",
    "# Define bin width and create bins\n",
    "binwidth = (max_cost - min_cost) / 30  # Adjust the divisor for desired number of bins\n",
    "bins = np.arange(min_cost, max_cost + binwidth, binwidth)\n",
    "\n",
    "# Plot the overlaid histograms with the same binwidth and range\n",
    "plt.figure(figsize=(10, 6))\n",
    "\n",
    "# Histogram for total cost using normal distribution costs\n",
    "sns.histplot(total_costs_df['Total_Cost_Sol1'], bins=bins, color='blue', label='Optimal Solution',\n",
    "             alpha=0.5, kde=False, stat='frequency')\n",
    "\n",
    "# Histogram for total cost using triangular distribution costs\n",
    "sns.histplot(total_costs_df['Total_Cost_Sol2'], bins=bins, color='orange', label='Second Best Solution',\n",
    "             alpha=0.5, kde=False, stat='frequency')\n",
    "\n",
    "# Add vertical lines at the means\n",
    "plt.axvline(total_costs_df['Total_Cost_Sol1'].mean(), color='blue', linestyle='dashed', linewidth=2)\n",
    "plt.axvline(total_costs_df['Total_Cost_Sol2'].mean(), color='orange', linestyle='dashed', linewidth=2)\n",
    "\n",
    "plt.title('Overlaid Histograms')\n",
    "plt.xlabel('Cle-Ind Cost')\n",
    "plt.ylabel('Frequency')\n",
    "plt.figtext(0.5, 0.01, 'Dashed lines show respective means', ha='center', fontsize=10)\n",
    "plt.legend()\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Looking at this histogram, how would you argue going forward with either solution?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Individual chance constraint\n",
    "\n",
    "Meeting widget demand in Nashville is now an extremely high priority. Historically, our data says that the demand here for the next sales period follows a [Poisson distribution](https://en.wikipedia.org/wiki/Poisson_distribution) that has a mean of 100 (in the original model, demand was 101). Fill in the code blanks below to make sure we meet the demand with probability 0.95."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "117.0"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from scipy.stats import poisson\n",
    "\n",
    "q95 = poisson.ppf(0.95, 100)\n",
    "q95"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use the demand constraint's attribute 'RHS' to quickly update the demand and use the code below to resolve and compare. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "This model has a total cost of 1798.3. The original model has a cost of 1748.42\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th>cost</th>\n",
       "      <th>shipment</th>\n",
       "      <th>shipment_Nash</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>production</th>\n",
       "      <th>distribution</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th rowspan=\"3\" valign=\"top\">Baltimore</th>\n",
       "      <th>Columbia</th>\n",
       "      <td>4.50</td>\n",
       "      <td>0.0</td>\n",
       "      <td>30.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Nashville</th>\n",
       "      <td>5.96</td>\n",
       "      <td>30.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Richmond</th>\n",
       "      <td>1.96</td>\n",
       "      <td>116.0</td>\n",
       "      <td>116.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"3\" valign=\"top\">Cleveland</th>\n",
       "      <th>Columbia</th>\n",
       "      <td>2.43</td>\n",
       "      <td>89.0</td>\n",
       "      <td>59.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Indianapolis</th>\n",
       "      <td>2.37</td>\n",
       "      <td>95.0</td>\n",
       "      <td>95.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Nashville</th>\n",
       "      <td>4.13</td>\n",
       "      <td>0.0</td>\n",
       "      <td>37.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Little Rock</th>\n",
       "      <th>St. Louis</th>\n",
       "      <td>2.92</td>\n",
       "      <td>140.0</td>\n",
       "      <td>140.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Birmingham</th>\n",
       "      <th>Nashville</th>\n",
       "      <td>1.53</td>\n",
       "      <td>71.0</td>\n",
       "      <td>80.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"2\" valign=\"top\">Charleston</th>\n",
       "      <th>Lexington</th>\n",
       "      <td>1.61</td>\n",
       "      <td>121.0</td>\n",
       "      <td>121.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>St. Louis</th>\n",
       "      <td>4.60</td>\n",
       "      <td>41.0</td>\n",
       "      <td>41.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                          cost  shipment  shipment_Nash\n",
       "production  distribution                               \n",
       "Baltimore   Columbia      4.50       0.0           30.0\n",
       "            Nashville     5.96      30.0            0.0\n",
       "            Richmond      1.96     116.0          116.0\n",
       "Cleveland   Columbia      2.43      89.0           59.0\n",
       "            Indianapolis  2.37      95.0           95.0\n",
       "            Nashville     4.13       0.0           37.0\n",
       "Little Rock St. Louis     2.92     140.0          140.0\n",
       "Birmingham  Nashville     1.53      71.0           80.0\n",
       "Charleston  Lexington     1.61     121.0          121.0\n",
       "            St. Louis     4.60      41.0           41.0"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "meet_demand['Nashville'].RHS = q95\n",
    "m.update()\n",
    "m.optimize()\n",
    "\n",
    "x_values_Nash = pd.Series(m.getAttr('X', x), name = 'shipment_Nash', index = transp_cost.index)\n",
    "\n",
    "m_Nash_ObjVal = m.ObjVal\n",
    "print(f\"This model has a total cost of {round(m_Nash_ObjVal,2)}. The original model has a cost of {round(m_ObjVal,2)}\")\n",
    "sol_Nash = pd.concat([transp_cost, x_values_Nash], axis=1)\n",
    "sol_Nash_comp = pd.concat([sol, sol_Nash.drop('cost', axis = 1)], axis=1)\n",
    "sol_Nash_comp[(sol_Nash_comp.shipment > 0) | (sol_Nash_comp.shipment_Nash > 0)]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Scenario Optimization\n",
    "\n",
    "Here we build on the example from the earlier presentation where we were considering three scenarios in portfolio optimization. We'll combine the concept of scenario optimization with \"maxi-min\" from Opti 201. We are considering investing in seven possible sectors under three different scenarios. Below is the input data which defines the sectors, scenarios, and the respective estimated return. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define sectors\n",
    "sectors = [\n",
    "    'Tech',\n",
    "    'TradRetailUtil',\n",
    "    'Renewable',\n",
    "    'FossilFuel',\n",
    "    'Auto',\n",
    "    'GenHealth',\n",
    "    'TradHealth'\n",
    "]\n",
    "\n",
    "# Define scenarios\n",
    "scenarios = [1, 2, 3]\n",
    "\n",
    "# Expected returns (%) for each sector under each scenario\n",
    "estimated_returns = {\n",
    "    ('Tech', 1): 12,\n",
    "    ('Tech', 2): 1,\n",
    "    ('Tech', 3): 3,\n",
    "    ('TradRetailUtil', 1): 3,\n",
    "    ('TradRetailUtil', 2): 4,\n",
    "    ('TradRetailUtil', 3): 3,\n",
    "    ('Renewable', 1): 5,\n",
    "    ('Renewable', 2): 15,\n",
    "    ('Renewable', 3): 4,\n",
    "    ('FossilFuel', 1): 3,\n",
    "    ('FossilFuel', 2): -5,\n",
    "    ('FossilFuel', 3): 4,\n",
    "    ('Auto', 1): 4,\n",
    "    ('Auto', 2): 2,\n",
    "    ('Auto', 3): 4,\n",
    "    ('GenHealth', 1): 1,\n",
    "    ('GenHealth', 2): 4,\n",
    "    ('GenHealth', 3): 10,\n",
    "    ('TradHealth', 1): 3,\n",
    "    ('TradHealth', 2): 2,\n",
    "    ('TradHealth', 3): -3\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create an optimization model where we decide the percentage of our budget to invest to **maximize the minimum** total estimated returns considering all scenarios. **HINT**: Define a single new auxiliary variable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimal objective value: 5.93\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Investment Allocation</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Tech</th>\n",
       "      <td>0.350</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>TradRetailUtil</th>\n",
       "      <td>0.000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Renewable</th>\n",
       "      <td>0.271</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>FossilFuel</th>\n",
       "      <td>0.000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Auto</th>\n",
       "      <td>0.000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>GenHealth</th>\n",
       "      <td>0.380</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>TradHealth</th>\n",
       "      <td>0.000</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                Investment Allocation\n",
       "Tech                            0.350\n",
       "TradRetailUtil                  0.000\n",
       "Renewable                       0.271\n",
       "FossilFuel                      0.000\n",
       "Auto                            0.000\n",
       "GenHealth                       0.380\n",
       "TradHealth                      0.000"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Create a new model\n",
    "m_scen = gp.Model(\"PortfolioOptimization\")\n",
    "m_scen.setParam('OutputFlag', 0)\n",
    "\n",
    "# Decision variables: Investment fractions in each sector\n",
    "x = m_scen.addVars(sectors, name=\"Investment\", lb=0, ub=1)\n",
    "\n",
    "# Auxiliary variable: Minimum portfolio return\n",
    "z = m_scen.addVar(name=\"MinReturn\")\n",
    "\n",
    "# Add investment allocation constraint\n",
    "m_scen.addConstr(x.sum() == 1, name=\"TotalInvestment\")\n",
    "\n",
    "# Add portfolio return constraints for each scenario\n",
    "portfolio_return = m_scen.addConstrs((gp.quicksum(estimated_returns[(i, s)] * x[i] for i in sectors) >= z for s in scenarios), name = 'ReturnConstraint_Scenario')\n",
    "\n",
    "# Set objective: Maximize the minimum portfolio return\n",
    "m_scen.setObjective(z, GRB.MAXIMIZE)\n",
    "\n",
    "# Optimize the model\n",
    "m_scen.optimize()\n",
    "\n",
    "print(f\"Optimal objective value: {round(m_scen.ObjVal,2)}\")\n",
    "\n",
    "# Get the results\n",
    "investment_allocation = pd.Series({i: round(x[i].X,3) for i in sectors}, name='Investment Allocation')\n",
    "investment_allocation.to_frame()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You now decide that your investments need to be a little more diversified. Add any necessary variables and constraints that will:\n",
    "- Make sure you invest in at least four sectors.\n",
    "- If an investment is made in a sector, then will need to be at least 10% of the budget."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Investment Allocation</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Tech</th>\n",
       "      <td>0.307</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>TradRetailUtil</th>\n",
       "      <td>0.000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Renewable</th>\n",
       "      <td>0.257</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>FossilFuel</th>\n",
       "      <td>0.000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Auto</th>\n",
       "      <td>0.100</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>GenHealth</th>\n",
       "      <td>0.336</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>TradHealth</th>\n",
       "      <td>0.000</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                Investment Allocation\n",
       "Tech                            0.307\n",
       "TradRetailUtil                  0.000\n",
       "Renewable                       0.257\n",
       "FossilFuel                      0.000\n",
       "Auto                            0.100\n",
       "GenHealth                       0.336\n",
       "TradHealth                      0.000"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Add new variables and/or constraints\n",
    "y = m_scen.addVars(sectors, vtype=GRB.BINARY,name='decision')\n",
    "\n",
    "l = 0.1\n",
    "m_scen.addConstrs((x[i] <= y[i] for i in sectors), name='xy_link1')\n",
    "m_scen.addConstrs((x[i] >= l*y[i] for i in sectors), name='xy_link2')\n",
    "m_scen.addConstr(y.sum() >= 4)\n",
    "# Optimize the model\n",
    "m_scen.optimize()\n",
    "\n",
    "# Get the results\n",
    "investment_allocation = pd.Series({i: round(x[i].X,3) for i in sectors}, name='Investment Allocation')\n",
    "investment_allocation.to_frame()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "gurobi_ml",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
